arXiv:1505.03320v2 [physics.data-an] 7 Jan 2016 


Integration of in situ Imaging and Chord Length Distribution 
Measurements for Estimation of Particle Size and Shape 


Okpeafoh S. Agimelen a ’*, Anna Jawor-Baczynska a , John McGintyg Christos 
Tachtatzis b , Jerzy Dziewierzg Ian Haley d , Jan Sefcik a , Anthony J. Mulholland 6 ’* 

a EPSRC Centre for Innovative Manufacturing in Continuous Manufacturing and Crystallisation, 
Department of Chemical and Process Engineering, University of Strathclyde, James Weir Building, 
7 5 Montrose Street, Glasgow, G1 1XJ, United Kingdom. 
b Centre for Intelligent Dynamic Communications, Department of Electronic and Electrical 
Engineering, University Of Strathclyde, Royal College Building, 204 George Street, Glasgow, Cl 

1XW, United Kingdom. 

c The Centre for Ultrasonic Engineering, Department of Electronic and Electrical Engineering, 
University Of Strathclyde, Royal College Building, 204 George Street, Glasgow, Cl 1XW, United 

Kingdom. 

d Mettler-Toledo Ltd., 64 Boston Road, Beaumont Leys Leicester, LE4 1AW, United Kingdom 
e Department of Mathematics and Statistics, University of Strathclyde, Livingstone Tower, 26 
Richmond Street, Glasgow Cl 1XH, United Kingdom 


Abstract 

Efficient processing of particulate products across various manufacturing steps re¬ 
quires that particles possess desired attributes such as size and shape. Controlling the 
particle production process to obtain required attributes will be greatly facilitated us¬ 
ing robust algorithms providing the size and shape information of the particles from 
in situ measurements. However, obtaining particle size and shape information in situ 
during manufacturing has been a big challenge. This is because the problem of estimat¬ 
ing particle size and shape (aspect ratio) from signals provided by in-line measuring 
tools is often ill posed, and therefore it calls for appropriate constraints to be imposed 
on the problem. One way to constrain uncertainty in estimation of particle size and 
shape from in-line measurements is to combine data from different measurements such 
as chord length distribution (CLD) and imaging. This paper presents two different 
methods for combining imaging and CLD data obtained with in-line tools in order to 
get reliable estimates of particle size distribution and aspect ratio, where the imaging 
data is used to constrain the search space for an aspect ratio from the CLD data. 

Keywords: Chord Length Distribution, Particle Size Distribution, Particle Shape, 
Focused Beam Reflectance Measurement, Imaging. 


* Corresponding authors 

Email addresses: okpeafoh.agimelen@strath.ac.uk (Okpeafoh S. Agimelen), 
anthony.mulholland@strath.ac.uk (Anthony J. Mulholland) 


Preprint submitted to Chemical Engineering Science 


January 8, 2016 



1. Introduction 


One of key steps in the manufacture of particulate products in the pharmaceuticals 
and fine chemicals industry is crystallisation, which is widely used for separation and 
purification of intermediates, fine chemicals and active pharmaceutical ingredients. The 
crystals come in different sizes and shapes. Subsequent steps in the manufacturing 
process, such as filtration, drying, blending and formulation of final products, require 
that the particle sizes and shapes lie within some desirable range. In order to provide 
monitoring and control of crystallisation processes it is necessary to develop techniques 
for estimating the shape and size distribution of particles in situ. There are a number 
of off line tools [1] that can be used to estimate the particle size distribution 1 (PSD) 
of crystals produced in a crystallisation process. However, of particular importance to 
the control of a crystallisation process are tools that can be used in situ. These tools 
should be suitable for estimating size and shape information of particles dispersed in 
a slurry without the need for sampling and/or dilution. Examples of such instruments 
are the focused beam reflectance measurement (FBRM), the three dimensional optical 
reflectance (3D-ORM) [2] and the particle vision and measurement (PVM) [3] sensors. 

In-line sensors such as FBRM and 3D-ORM measure a chord length distribution 
(CLD) 2 which is related to the size and shape of the particles in a slurry. It has been 
a long standing challenge to be able to deduce the actual PSD and particle shape from 
experimental CLD data. In order to do this, an inverse problem needs to be solved. This 
is usually achieved by suitably discretising both CLD (which is already measured as 
a discrete distribution) and PSD and then constructing an appropriate transformation 
matrix relating these two distributions [4-6]. The transformation matrix depends on 
the choice of size bins used to discretise the two distributions and the corresponding 
size ranges as well as the shape of particles. The transformation matrix is usually 
not known in advance and needs to be estimated along with the corresponding PSD 
(discretised) so that the convolution of the transformation matrix with the PSD yields 
a CLD which agrees with the experimentally measured CLD. However, this problem 
is ill posed. There are a number of different transformation matrices and PSDs whose 
convolutions give rise to the same CLD. Hence the challenge is how to estimate a 
combination of transformation matrix and PSD whose convolution will agree with an 
experimentally measured CLD for a given slurry as well as the PSD estimated being 
physically reasonable and representative of the particles in the slurry. 

The approach which was used in previous works [7-10] when estimating the trans¬ 
formation matrix was to assume the same shape (quantified by a metric referred to as 
aspect ratio) for all the particles in the slurry, and then use a previously estimated 3 


1 The term particle size distribution is broadly used here to refer both to continuous analytical 
probability density functions for particle sizes and discretised probability histograms of the particle 
sizes. 

2 Similar to the case of PSD, the term chord length distribution is used to cover both continuous 
analytical functions and discretised probability histograms. 

3 The approach was to estimate the range of particle sizes in a sample by techniques such as sieving, 
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range of particle sizes in the slurry to construct the transformation matrix. This ap¬ 
proach is not suitable for monitoring a crystallisation process where nucleation and/or 
growth of particles was present as neither the range of particle sizes nor their aspect 
ratio would be known in advance. A technique which was suitable for estimating the 
range of particle sizes in a slurry in situ was presented in our previous work [6]. How¬ 
ever, like in other previous works, the transformation matrix was constructed with a 
single aspect ratio for all particles. This leaves open a possibility that the transforma¬ 
tion matrix is constructed with inappropriate aspect ratio or that there is a wider range 
of aspect ratios present for particles of same or different sizes. It was demonstrated in 
our previous work [6] that it was still possible to calculate different CLDs that all had 
a very good agreement with an experimentally measured CLD even though some of the 
transformation matrices were constructed at aspect ratios that were far from the shape 
of the particles described. However, it was also shown that as the aspect ratio deviated 
further from the true shape of the particles, then the corresponding PSD became in¬ 
creasingly noisy. This situation led to the introduction of a penalty function in order to 
eliminate unrealistic aspect ratios. However, when there is a wide variation of aspect 
ratios of the particles in the slurry, there is a need to introduce further constraints on 
the aspect ratio to reduce the search space and regularize the inverse problem. One 
way to do this is to get estimates of aspect ratio (within some reasonable bounds) using 
imaging, and then use this information to constrain the search for a representative as¬ 
pect ratio. However, the imaging needs to be done in situ in order to develop techniques 
for estimation of PSD and particle shape which could be used for real time monitoring 
and control of particle production processes. 

While it would be desirable to get good estimates of both PSD and particle aspect 
ratio using in situ imaging alone, this is currently not the case. The currently available 
in-line imaging tools (for example, the PVM used in this work) produce 2D projection 
images. Furthermore, the objects in the images may be partially or completely out of 
focus, parts of imaged object may cross the image frame or objects may overlap each 
other * * 4 . Although advanced measurement equipment have been developed which can be 
used to capture 3D images of particles in a slurry and make good estimates of PSD and 
shape of particles, it requires sampling and dilution flow loops 5 to allow capturing 3D 
images of individual particles in a flow-through cell [11-14], Therefore this approach 
may not be generally applicable for in-line monitoring of particle manufacturing pro¬ 
cesses. Hence the current situation is that PSD cannot be estimated to a good degree 
of accuracy using routinely available in-line imaging tools. To overcome this challenge, 
we propose to combine in-line CLD measurements with imaging data to provide more 
reliable estimate of quantitative particle attributes. 


laser diffraction, microscopy or use information supplied by the manufacturer before suspending the 

particles. 

4 The issue of objects overlapping each other would not be a problem if an appropriate image 
processing algorithm which can resolve the objects is used. 

5 The dilution is necessary to avoid instances of overlapping particles in images. 
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In this paper, we present two different methods for combining imaging data with 
CLD data for particle size and shape estimation. The first method presented here 
calculates an estimate of the mean aspect ratio of all the particles in the slurry and 
then uses this information to constrain the search space for size and shape estimation 
from the CLD data. In the second method, a distribution of aspect ratios for each 
particle size is used for the PSD estimation. The distribution of aspect ratios is based 
on the data from the captured images. 

2. Experiment 

To demonstrate the technique for estimating particle size and shape information 
using a combination of the CLD and imaging data, experiments were performed in 
slurries containing particles of different shapes. The materials and procedure are de¬ 
scribed below. 

2.1. Materials 

Three different samples were used for the measurements. Sample 1 consisted of 
polystyrene (PS) microspheres purchased from EPRUI nanoparticles and Microspheres 
Co. Ltd. with batch number 2012-5-7, and 0.2g of the PS microspheres were dispersed 
in lOOg of isopropanol (IPA) purchased from VWR (20842.323) giving a concentra¬ 
tion of 0.2% by weight. Sample 2 consisted of cellobiose octaacetate (COA) particles 
obtained from GSK. The particles were dispersed in methanol (purchased from VWR 
(20847.307)) with the same concentration as in sample 1. Sample 3 consisted of glycine 
(Glycine) crystals obtained by cooling crystallisation from an aqueous solution. The 
solution with glycine concentration of 340mg/ml was prepared using glycine (purchased 
from Sigma-Aldrich (G8898, > 99% TLC)) and deionized water (from an in-house Mil- 
lipore Water System (18Mf2/cm)). The solution was cooled from a temperature of 
90°C to a temperature of 43°C at a rate of 3°C/min. During this process the glycine 
crystallised out of solution until an equilibrium particle size distribution was reached. 
The crystallisation of glycine was monitored with the FBRM probe which showed an 
initial increase (in time) of chord lengths before eventually reaching a steady state. 

2.2. Experimental Setup 

The suspension of particles for all samples was made in the Mettler Toledo EasyMax 
102 system. The EasyMax system consists of a cylindrical jacketed vessel of volume 
100ml with different stirrer and blade options. An anchored overhead stirrer with 
pitched (45° pitch angle) blades was used in all the experiments in this work. The 
stirrer shaft and probes were inserted into the slurries through ports located at the top 
of the set up. The stirring speed was set at 400 rpm in all experiments. 

The CLD measurements were made with a Mettler Toledo FBRM G400 probe. The 
images of the particles were captured with a Mettler Toledo PVM V819 probe during 
the period of CLD measurement. The FBRM sensor consists of a system of lenses which 
focus a laser beam onto a spot near the probe window in the slurry. The laser spot 
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Figure 1: The CLD (measured with the FBRM G400 sensor) for the (a) PS, (b) COA and (C) Glycine 
samples. 


moves in a circular trajectory and the back scattered light is detected. The chord length 
is then calculated as the speed of the laser spot multiplied by the duration of the back 
scattered light as the laser traverses a particle. The FBRM sensor records the lengths 
of the chords for a pre-set duration after which the CLD is reported [2, 6, 15-17]. 

The PVM is an in situ microscope which consists of eight laser sources enclosed in 
a cylindrical tube. The six forward lasers and two back lasers (achieved by reflecting 
two lasers off a Teflon cap at the probe window) illuminate the particles in the slurry. 
The back scattered light is detected by a CCD element from which grayscale images are 
constructed. The image frame of the CCD array consists of 1360x1024 pixels with a 
pixel size of 0.8/irn. The PVM V819 sensor has a maximum acquisition rate of 5 images 
per second, although lower rates of acquisition could be set depending on requirements. 
The depth of the focal zone is restricted to about 50/irn so that all objects that are in 
focus result in images that have identical magnification levels. Each of the lasers can 
be switched on or off so that different degrees of illumination can be achieved. 

3. Experimental Data 

The CLD data measured with the FBRM sensor for PS is shown in Fig. 1(a). Images 
were captured with the PVM sensor during the CLD measurement. A representative 
image from the PS sample is shown in Fig. 2(a). 

Similarly, the CLD data recorded for COA is shown in Fig. 1(b) while that of 
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Figure 2: Representative images (obtained with the PVM V819 sensor) for the (a) PS, (b) COA and 
(c) Glycine samples. The images have the same width of 1088/nn and height of 819^m. 


Glycine is shown in Fig. 1(c). Representative images for COA and Glycine are shown 
in Figs. 2(b) and 2(c) respectively. A total of 600 images were acquired each for PS 
and COA, and 121 images for Glycine. 

4. Image Analysis 

As mentioned in the introductory section, an estimate of the PSD and particle shape 
can be made from images alone without the need to include CLD data. However, due to 
the reasons discussed in the introductory section, this approach is not always convenient. 
The techniques presented here utilise images obtained with an in-line measuring tool. 
However, due to the limitations in the images (as discussed in the introductory section), 
it is necessary to combine the imaging data with CLD to obtain reasonable aspect ratio 
and/or size estimates. 

The images captured with the PVM are processed in order to detect the objects 
contained in them and hence obtain information about the shape of the particles in 
the slurry. However, the image processing algorithm used in this work does not have 
features to resolve an object completely when it does not lie entirely in focal plane of 
the PVM sensor. Also, it does not have functionalities to resolve overlapping objects in 
images. For this reason the samples used in this work were deliberately prepared dilute 6 


6 Low slurry densities have been used here for the purpose of methods development. Future work will 
involve the investigation of the applicability of the methods developed here at higher slurry densities 
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to reduce instances of overlapping objects. The parameters' of the image processing 
algorithm were tuned to reject most of the particles that were not entirely in the focal 
plane of the PVM sensor. However, the imaging data still has some degree of inac¬ 
curacy as seen in the error bars of the data (see subsection 4.6 of the supplementary 
information). This limitation not withstanding, the data from the images was suffi¬ 
ciently accurate to demonstrate the techniques developed in this work. Furthermore, 
images which do not contain objects that are contained completely within the image 
frame were also discarded. This situation of having to discard some images reduces the 
number of data sets that can be gathered from the images. However, it can be shown 
(see section 2 of the supplementary information) that with a sample size (number of 
objects) of about 500 the error incurred in estimating the aspect ratio is reduced to a 
reasonable extent. However, for a more robust estimate of the PSD using imaging data, 
a larger number of objects will be required. The number * * * * 8 of detected objects used in 
this work is just sufficient to demonstrate the methods developed here. The issue of 
objects not completely in focus can be dealt with if additional functionalities are added 
to the image processing algorithm, but this is beyond the scope of this work. The key 
steps for detecting objects in the images captured by the PVM sensor are summarised 
in subsection 4.1. 

4-1. Object Detection 

The raw grayscale images from the PVM sensor are passed through a median filter 
to remove speck noise from the image background which is homogeneous. At this stage 
objects on the boundary of the image frame are removed. Any object with surface area 
below 900 pixels 9 is considered noise and excluded from processing. Finally, a closing 
operation with a disk structural element is used to join broken edges. The resulting 
blob properties such as area, centroid, eccentricity, convex area, and major and minor 
axes can be obtained. These steps are summarised in Fig. 3. The tunable parameters of 
the image processing algorithm are summarised in subsection 4.5 of the supplementary 
information. 

4-2. Characterising Particle Shape 

The following procedure is used to obtain a shape descriptor, which is then used to 
characterise the shape of each particle. Each boundary pixel j of object i has coordinates 


(with a more advanced image processing algorithm) using suspensions of particles of known PSD, 

then the degree of deviation of the results from the known PSD can be quantified at different slurry 

densities. 

'The parameters of the algorithm need to be tuned for different samples due to variation of contrast. 

8 A total of 1393 objects were detected for PS, 1810 for COA and 526 for Glycine. 

9 The surface area of 900 pixels represents length dimensions of approximately 30 x 30 pixels 
(assuming a square geometry). This implies that objects that are smaller than approximately 24/rm 
are rejected by the image processing algorithm. The consequence is that there is no estimate of aspect 
ratio for these small objects. However, the particles used in this work have sizes mostly in the range 
of 100/rm so that the effects of this are minimal. 


7 




Figure 3: A schematic for the key steps of the image processing algorithm used in this work. 




Figure 4: A sketch of the particle characterisation procedure. The boundary pixels before rotation are 
shown in (a). The line joining the farthest boundary pixel e tJ * to the centroid makes an angle p> 
with the horizontal axis, (b) A rotation is performed so that the farthest pixel from the centroid lies 
on the horizontal axis. The indices of the remaining pixels are assigned relative to this pixel. 


e ij = [e lJt xi &ij,y], where the Cartesian coordinates system has been used. The centroid 
of object i has coordinates e,. The distance d,j of the centroid of object i to the pixel 
j is given as 

dij 11 e^j | . (1) 

Then for each object i, locate the pixel (with index j*) whose distance from the centroid 
d^* is largest as 

d^* = max {dij}. (2) 

j 

The angle between the line joining this farthest pixel to the centroid and the horizontal 
axis is ip as shown in Fig. 4(a). Then, for each object i, transform the coordinates of 
each pixel by performing a rotation through the angle ip so that the line joining the 
centroid to the farthest pixel becomes parallel to the horizontal axis. This operation 
transforms the coordinates of each pixel j to e t] and the centroid to §j. The line joining 
the pixel j to the centroid now makes an angle 6 with the horizontal axis as shown 
in Fig. 4(b). Due to the rotation, the farthest pixel from the centroid corresponds to 
9 = 0. The distance of each pixel j from the centroid is the same of course. 
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Figure 5: The shape descriptor (d in Eq. (4)) as a function of the angle (9) for the (a) PS, (b) COA 
and (c) Glycine samples. 


Since the sample rate with respect to 9 is not uniform, then the pixel distances 
were resampled with N p uniformly spaced 6 values constructed as 9 P = pA9, A9 = 
2i t/N p , p = 1,2 ,..., N p . This allows the vector of all pixels distances (from the 
centroid) for object i to be written as 

d* [dy i. di2j ■ ■ ■, djTvp]- (3) 


Finally, the average vector d of all pixel distances for all objects detected by the image 
processing algorithm can be calculated as 


d 


1 

N 0 bj 


N 0 bj 




( 4 ) 


where N a bj is the number of objects detected from all the images analysed. 

A plot of d versus angle (6) can be made as shown in Figs. 5(a) to 5(c). For near 
spherical particles, the shape descriptor d is nearly constant as in the case of PS in 
Fig. 5(a). However, for elongated particles, the shape descriptor d has two minima 
and maxima as in the cases of CoA and Glycine in Figs. 5(b) and 5(c). The shape 
descriptor shown in 5(a) to 5(c) is similar to the type described in [13]. 

In an ideal situation, the shape descriptor for spherical particles will be constant 
at a value representing the radius of the spherical particles. However, since the PS 
particles are not perfectly spherical there is slight variation in the dimensions so that 
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an average aspect ratio r (the ratio of the minor to the major dimension) can be 
estimated. Similarly, the maxima in the shape descriptors for the CoA and Glycine 
particles (in Figs. 5(b) and 5(c)) represent the major dimension of the particles while 
the minima represent the minor dimension of the particles. 

For the case of the PS particles the mean aspect ratio r is estimated from the 
minimum dimension C m i n and maximum dimension C max (see Fig. 5(a)) as 


Cr 


r = 


Cr 


(5) 


In the case of elongated particles, the maximum dimension (average length of particles) 
is given as C max = Ci max + C 2 max (shown in Figs. 5(b) and 5(c)) and the minimum 
dimension (average width of particles) is given as C nun = C\ min + £ 2 mm- So that the 
average aspect ratio can then be estimated using Eq. (5). 

However, the aspect ratios for individual particles will be different from r. The 
aspect ratio for each particle is estimated from the shape descriptor corresponding 
to that particle. Once the aspect ratios of individual particles are estimated, then a 
scatter plot of aspect ratio versus particle length can be made as shown in Figs. 6(a) 
to 6(c). The shape descriptors for individual particles are not always smooth as in 
the cases shown in Figs. 5(a) to 5(c). They contain different degrees of variation due 
to imperfections in the particles and images (see subsection 4.6 of the supplementary 
information for details). 

Furthermore, the aspect ratios estimated from the shape descriptor for individual 
particles contain some artefacts (as can be seen in Fig. 6) due to a number of factors. 
These factors include deformations in the particles (that is, particles whose shapes 
deviate from the majority of particles in the slurry), impurity objects in the monitored 
slurry and particles not completely in focus 10 . For example, aspect ratios as low as 
about 0.2 in the case of the spherical particles in Fig. 6(a) or the minor peak at aspect 
ratio r & 0.5 in Fig. 6(d) are artefacts. 

Figures 6(b) and 6(c) suggest the presence of particles of sizes up to about 500pm 
in the COA and Glycine samples respectively. However, the number of data points (in 
Figs. 6(b) and 6(c)) corresponding to these large particles may not be representative 
of the actual number of these large particles in the slurry. This is because the image 
processing algorithm has been designed to remove objects making contact with the 
image frame, and larger particles have a higher probability of making contact with 
the image frame. In the situation where the PSD is to be estimated from image data 
alone, then this probability will need to be taken into account [18]. However, since the 
objective here is just to estimate the aspect ratio of particles of different sizes, then this 
is not a crucial issue. 


10 The image processing algorithm parameters are tuned to remove particles that are out of focus. 
However, when the particle is only partially in focus, the image processing algorithm detects only part 
of the particle and this leads to an error in the estimated aspect ratio for that particle. 


10 





Particle length(^m) Particle length(pm) 


Particle length (pm) 



Figure 6: Scatter plots of the aspect ratio versus particle length for individual objects for (a) PS, 
(b) COA and (c) Glycine. The histograms of the aspect ratios are shown in (d) PS, (e) COA and 
(f) Glycine. 


The histogram (in Fig. 6(d)) 11 for the spherical particles has a dominant mode close 
to the aspect ratio of 0.85 which is close to the average aspect ratio r ~ 0.8 estimated 
from the average shape descriptor (in Fig. 5(a)) for this slurry. Similarly the modes of 
the histograms (in Figs. 6(e) and 6(f)) for COA and Glycine occur close to aspect ratios 
of 0.2 and 0.4 respectively. These values are close to the average aspect ratios r ~ 0.2 
(for COA) and r « 0.4 (for Glycine) obtained from their respective shape descriptors 
in Figs. 5(b) and 5(c). 

5. Modelling Chord Length Distribution 

The sizes of particles in a slurry can be represented by the equivalent spherical 
diameter as was done in [6]. However, a characteristic length L could also be used, 
which can be chosen as the distance between the two extreme points in the particles’ 
geometry. Since the estimated sizes from the images is L, then this metric is used here 
for consistency with the image data. Once the metric for particle sizes has been chosen, 
then the PSD can then be expressed in terms of the chosen particle size metric. The 


11 The uniform bin widths of the histograms in Figs. 6(d) to 6(f) were estimated using the Freedman- 
Diaconis rule [19]. The number of bins were then estimated as 18 bins for PS, 28 bins for COA and 
13 bins for Glycine. 
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PSD X is related to the CLD C by means of a convolution function [2, 6, 20] and the 
relationship can be written in matrix form as [6] 12 

C(s) = A(s, L)X(L), (6) 

where s is the chord length, and X is the length weighted PSD [6] given as 

Xi = LiXi, i = 1,2,3, N. (7) 

The PSD X t (which is actually a histogram) consists of N bins. The characteristic size 
Li of particle size bin i is the geometric mean of sizes and L i+ i as L { = -y/LjLj +1 . 
The bin boundaries are calculated as 


Li — L min uj l 1 , i — 1, 2,..., TV + 1 


( 8 ) 


where 


oj = 


N 


1 


( 9 ) 


where L m i n is the left boundary of the first particle size bin and L max is the right 
boundary of the last particle size bin. 

In previous works [7-9, 21-23] the values of L min and L max were estimated from 
suitable measurements. However, the technique of estimating L min and L max directly 
from the bin boundaries of the CLD histogram using a moving window technique has 
been demonstrated to yield more accurate results [6]. This window technique is more 
suitable for estimating the sizes of particles in-line in a process where particle size 
information is obtained from the CLD [6]. 

The length weighting applied to the PSD X in Eq. (7) is necessary because the CLD 
for a population of particles is biased towards particles of larger sizes [6, 20, 24, 25]. The 
forward problem in Eq. (6) is implemented by considering a chord length histogram Cj 
of M bins where the characteristic chord length Sj of bin j is the geometric mean of 
the chord lengths of s 3 and s J+ \ as outlined in [5, 6]. 

If the PSD for a population of particles is known, then the CLD can be calculated 
using Eq. (6). However, in practical situations, the particle size histogram X,; is not 
known in advance resulting in the inverse problem of calculating an unknown PSD X t 
from a known CLD Cj. For this reason, the forward problem in Eq. (6) is reformulated 
as 

C(s) = A(s, L)X(L), (10) 

where the matrix A is obtained from matrix A by multiplying each column of A by 
the corresponding particle length as described in [6]. 


12 Note that the symbol L was used to represent the length of a chord in [6]. However, the symbol 
L is used to represent the characteristic length and s the length of a chord in this work. The CLD 
C and PSD X in Eq. (2) have been discretised. As such they are not continuous probability density 
functions and the term distribution is used in this work for simplicity. 
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Each column i of matrix A. ]r is calculated from the CLD of a single particle (single 
particle CLD) of length Li and given aspect ratio ty. In the current work, the single 
particle CLD used in constructing the columns of matrix A are obtained from the 
analytical Li and Wilkinson (LW) model [5] 13 . 

The process of calculating the single particle CLD involves computing the relative 
likelihood of obtaining a chord of length s from a particle of a given length and aspect 
ratio [5]. The LW model gives a probability density function (PDF) which can be 
used in making this calculation for ellipsoidal shaped particles 11 . The PDF is derived 
from an ellipsoid of semi major axis length a, semi minor axis length b and aspect ratio 
r = b/a [5]. The LW model gives the probability PzMji s j+i) °f obtaining a chord whose 
length lies between Sj and s^+i from an ellipsoid of characteristic length Lj = 2a* (see 
section 1 of the supplementary information and [5, 6] for the mathematical expression 
for pj j .( y Sj,Sj + 1 )). Once the probabilities are calculated, then for each row j of matrix 
A the columns are constructed as 

A j = [PL 1 ( s j, s j+i),PL 2 ( s j, s j+ 1), • • • ,PlS s ^ s 3+ 1), • • -,PL N { s j> s j+ 1)] • (H) 

6. Incorporating aspect ratio from images 

As stated in section 5 the calculation of a column of matrix A requires the charac¬ 
teristic size of the particle size bin corresponding to that column, as well as the aspect 
ratio r of the particle of that characteristic size. In the previous work [6], all parti¬ 
cles were assumed to have the same aspect ratio, and its value was estimated using an 
algorithm based solely on CLD data. The single aspect ratio approach will be used 
in subsection 6.1, where the single aspect ratio value is estimated from imaging data. 
This approach is most suitable for the case of spherical particles where the aspect ra¬ 
tios of the individual particles are tightly packed around some mean value. However, 
for the case of particles where there is a wider spread of aspect ratios, a variation of 
aspect ratios for different particle sizes can also be used. The corresponding technique 
is outlined in subsection 6.2. 

6.1. Method 1: Population CLD with a single aspect ratio 

When the aspect ratios are tightly packed around some mean value as in the case 
of spherical (or near spherical) particles in Figs. 6(a) and 6(d), it may be desirable 


13 Even though the model by Vaccaro et al. [25] gave estimates of particle aspect ratios that were 
closer to the estimates from images in [6], the LW model is used in this work for the CLD calculation. 
The reason is that the Vaccaro model is restricted to small values of aspect ratios r < 0.4, whereas the 
image data in Fig. 3 cover aspect ratios of r ~ 0.1 torw 0.9. The LW model covers the entire range 
from r = 0 to r = 1. 

14 The shapes of the COA and Glycine particles in this work have been approximated as ellipsoids. 
However, this is only an approximation as Fig. 2(c) clearly shows that the Glycine particles are 
faceted. The use of ellipsoids to represent faceted objects introduces some discrepancy between the 
single particle CLDs of both objects (see section 8 of the supplementary information for details). 
However, the ellipsoid approximation used here is sufficient to illustrate the methods presented here. 
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to use a mean value for all particles. Although, the mean aspect ratio estimated from 
the images provides the best estimate based on available data, there is an uncertainty 
due to sampling limitations and various artefacts discussed above. This necessitates a 
search (within a suitable confidence interval) around the mean aspect ratio estimated 
from the images for an aspect ratio which best matches the experimentally measured 
CLD. This results in the search space being narrowed down leading to significantly 
lower computation time and less uncertainty in subsequent calculations. 

In this section, the transformation matrix A in Eq. (10) is constructed with a single 
aspect ratio for all its columns. The aspect ratio r is chosen from a range given as 

r e [f — A fa r , f + A fcr r \, (12) 

where ay is the standard deviation of all aspect ratios estimated from images and J\f is 
the number of standard deviations chosen. The purpose of Eq. (12) is to constrain the 
search space for the inverse problem. The details of the inverse problem will be given in 
section 7. Hence, Method 1 corresponds to the case described in [6] where each column 
of matrix A consists of the single particle CLD of a particle of length L and aspect 
ratio r, where the search space for r is reduced by means of Eq. (12). 


6.2. Method 2: Population CLD with multiple aspect ratios 


The Method 1 presented in the previous subsection assigns the same aspect ratio to 
all the particles in the slurry. This method is capable of getting reasonable estimates 
of the PSD. However, to take the variation of particle shape into account, a second 
method is presented here in which multiple aspect ratios are assigned to particles of the 
same size. This method is particularly relevant for particles whose shape is needle-like 
or near rectangular as illustrated on the left of Fig. 7. The second method is outlined 
below. 

In Method 2, the particles in bin i are assigned the same characteristic size Lj but 
different aspect ratios. The aspect ratios assigned to the particles run from r" im to 
r m ax ag iu us t ra |;ed in Fig. 7. The procedure is to divide the particles in bin i into N r 
subgroups (N r = 50 in this case, see subsection 4.1 of the supplementary information 
for details), and the particles in subgroup k are assigned the same aspect ratio rf. The 
aspect ratios of the subgroups are uniformly spaced so that the aspect ratios are given 
as 


r 


e { r i = r fv‘ 


+ Ar,.. 


,fc N r max \ 

i i ■■■■>' i i J > 


(13) 


where 

r max _ T min 

Ar = ^-L-. (14) 

However, the aspect ratios r\ of the different subgroups need to be weighted by 
different amounts as the numbers of particles with given aspect ratios in the images 
may not necessarily be uniform from r™ ,m to r™ a3; . The data on the right of Fig. 7 
suggests an approximately uniform distribution of the aspect ratios from r™” to r™ ax , 
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Figure 7: Left: A schematic to illustrate the assignment of different aspect ratios to particles of the 
same length. Right: A scatter plot to illustrate the maximum aspect ratio r ™ a:r and minimum aspect 
ratio for a particular bin i for calculating the columns (of slice i) of the 3 D matrix in Eq. (15). 


hence the aspect ratios were assigned to each of the subgroup k with equal weight in 
this case. 

Once aspect ratios have been assigned to different particle size bins 15 , then a 3D 
probability array Ajki is constructed. Each slice i of the array corresponds to a par¬ 
ticle of size Li, and each column k of slice i contains the probabilities p T k(sj,Sj + 1 ) of 

obtaining chords whose lengths lie between Sj and Sj + i from a particle of size Lj and 
aspect ratio rf. Hence the 3D array consists of M rows, N r columns and N slices. The 
transformation matrix A,-, in Eq. (10) is then obtained from the 3D array by averaging 
over the slices as 

^ N r 

Ah = Ajki- (15) 

r k =1 

This simple averaging is carried out since the aspect ratios rf are assigned to the 
subgroups k with equal weights. This is the simplest way to construct the matrix 
from the 3D array Ajki, an d this approach is supported by the data from the images as 
shown on the right of Fig. 7. It is possible to introduce a probability distribution for the 
aspect ratios assigned to the subgroups, but this simple approach has been used here 
for the purpose of illustrating the technique. Once the matrix Aji has been constructed, 
then the forward problem in Eq. (10) can be solved for a given PSD X t . 


15 Details of the technique for combining the aspect ratios with the windowing technique developed 
in [6] can be found in section 5 of the supplementary information. 
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7. PSD Estimation 


As mentioned in section 5, the problem encountered in practical situations is the 
estimation of the PSD X t corresponding to an experimentally measured CLD C*. This 
is the inverse problem to the forward problem given in Eq. (10). One of the key steps in 
the process of the PSD estimation is to determine the transformation matrix A in Eq. 
(10) as accurately as possible. The level of accuracy of the matrix A depends on the 
values of L min and L max as well as the aspect ratio(s) used in calculating its columns. 

To determine the best possible values of L min and L max , the forward problem in Eq. 
( 10 ) is rewritten as 

C = AX + e, (16) 

where e is an additive error between the calculated and experimentally measured CLD. 
Then for given values of L min , L max , values of the fitting parameter 7 are found 16 which 
minimises the objective function /1 given as 


M 


/. = £ 


j = 1 


N 


c; - £ Tv 


i= 1 


(17) 


where 

Xi = e 7i ,i = 1,2,3, (18) 

A trial solution of 7 * = 0 was used in the calculation of the vector X j from Eq. (17). 
Once the solution vector X t is obtained, then it is used to solve the forward problem 
in Eq. (10) to obtain a calculated CLD Cj, j = 1, 2,..., M, where M is the number of 
bins in the CLD histogram 11 . The procedure is repeated until the optimum values of 
L min and L max are found for which there is the best match between the calculated and 
experimentally measured CLD. 

The objective function f\ given in Eq. (17) is suitable for estimating the optimum 
values of L min and L max whether the same aspect ratio is assigned to all particles or a 
distribution of aspect ratios is assigned to particles of the same size. In the case where 
a single aspect ratio is assigned to all particles, the objective function /j is not suitable 
for picking out the best aspect ratio within the confidence interval of aspect ratios. This 


16 The Levenberg-Marquardt (LM) algorithm as implemented in Matlab was used in this work to 
solve the optimisation problem here. The PSD Xi is estimated by means of the parameter yj. An 
initial value of 7 j is passed on to the LM algorithm which then searches for the optimum value of 7 ., 
to fit the given CLD. Since the PSD Xi is defined as an exponential function in Eq. 18, then the 
parameter 7 , can take values in the interval (— 00 , + 00 ) and still give Xi > 0. This implies that the 
non negativity requirement on the PSD is maintained by the formulation of Xi in Eq. 18. Therefore 
the LM algorithm was run without the use of lower or upper bounds as the parameter 7 i is defined in 
(— 00 , + 00 ). 

17 The value of M = 100 was used in all the calculations here to mimic the number of bins set in the 
FBRM G400 sensor. The values of N = 70 and N = 50 were used in Methods 1 and 2 respectively. 
See subsection 4.2 of the supplementary information for more details on the choice of the values of N 
for the two methods. 
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task is accomplished with another objective function f 2 (to be introduced in subsection 
7.1). The calculation for estimating the best aspect ratio (within the confidence interval 
of aspect ratios) using the objective function f 2 is carried out with the values of L min 
and L max estimated with the objective function f\. However, in the other case where 
a distribution of aspect ratios is assigned to particles of the same size, the objective 
function f 2 is not used as the problem of determining the best aspect ratio has been 
removed. 

Once the optimum values of L min and L max and (in the case of Method 1) the best 
aspect ratio have been estimated, then the PSDs (both number and volume based) 
can be calculated. However, these PSDs may not be reasonably smooth, showing non¬ 
physical oscillations as is often the case when solving ill-posed problems. In such cases, 
a third objective function / 3 (to be introduced in subsection 7.3) is used to calculate 
smooth PSDs. The calculation is carried out using the optimum values of L min and L max 
obtained with the objective function f\ and in the case of Method 1, the best aspect 
ratio obtained with the objective function f 2 . The calculation of the smooth PSDs with 
the objective function / 3 is done using suitable criteria described in subsection 4.4 of 
the supplementary information. 

As stated above, the objective function f\ is used to obtain the optimum values of 
L min and L max for both Methods 1 and 2. A given pair of L min and L max are said to 
be optimum when the corresponding calculated CLD C has the best match with the 
experimentally measured CLD C*. The level of agreement between the calculated and 
experimentally measured CLD is assessed by computing the L 2 norm 

uc*-cii = yy us) 

The values of L m i n and L max for which the L 2 norm in Eq. (19) reaches a minimum 
are chosen as the optimum values. 


7. 1. PSD estimation for Method 1 

In Method 1 the search for the optimum values of L min and L max using the objective 
function f\ is done at each aspect ratio within the confidence interval in Eq. (12). 
However, for particles of a given shape, the L 2 norm in Eq. (19) initially decreases with 
increasing aspect ratio and then becomes level (see subsection 4.3 of the supplementary 
information for details). This leads to non-uniqueness in determining the optimum 
value of r [6]. This problem of non-uniqueness is removed by using a modified objective 
function / 2 which contains a penalty term to control the size of the calculated PSD 
vector as 
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N 
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f2 = 


3 = 1 


C* 


i=l 


AjiXi 


+ v XI v 2 ■ 


( 20 ) 


2=1 


where the parameter Ai sets the level of imposed penalty. The value of Ai is chosen 
by comparing the magnitudes of the terms in Eq. (20) (see subsection 4.3 of the 
supplementary information). The aspect ratio at which the objective function / 2 reaches 
a minimum is then chosen as the optimum. 
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The solution vector (which is a number based PSD) obtained from Eq. (20) is not 
necessarily smooth as the penalty imposed on the solution vector only restricts its mag¬ 
nitude. With this penalty function, the LM algorithm could settle on a solution vector 
that contains some local fluctuations but whose value of f 2 is slightly less than a nearby 
solution that is smooth. For this reason, a new objective function / 3 (see subsection 
7.3 for details) which contains a penalty term to control the second derivative (to im¬ 
prove the smoothness of the solution vector) of the solution vector is used to estimate a 
number based PSD whose corresponding CLD is compared with the experimental data. 


1.2. PSD estimation for Method 2 

In Method 2 (described in subsection 6.2) particles of different characteristic sizes 
Li are assigned a range of aspect ratios as outlined in subsection 6.2. This eliminates 
the need to search for the best global aspect ratio which is the situation in Method 1. 
Hence in Method 2, it is only necessary to search for the best values of L min and L max 
using the objective function fi. Once the optimum transformation matrix is obtained 
using the objective function /i, then the corresponding smoothed solution is obtained 
with the objective function / 3 (given in Eq. (21)). 


7.3. Volume based PSD 

It is often necessary to recast the PSD X t (which is number based) in Eq. (17) 
as a volume based PSD since most instruments for measuring PSD give the data in 
terms of a volume based PSD. A new technique which allows suitable penalties to be 
imposed on the calculated volume based PSD X 1 ' was introduced in [6]. In the current 
work, a smoothing penalty (referred to in subsection 7.1, see also [26]) is imposed. This 
is because the estimated volume based PSD X’’ may contain significant non-physical 
oscillations even though the corresponding number based PSD only contains none or 
minor fluctuations. The objective function / 3 (see Eq. (21)) used to impose smoothness 
on the volume based PSD can also be used to obtain a smooth number based PSD. 
Hence the function / 3 is given in Eq. (21) in terms of generic quantities depending on 
whether the number based or volume based PSD is being computed. 

The function / 3 is given as 
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( 21 ) 


In the case of a number based PSD, the CLD Cj~ = C* (the experimentally measured 
CLD), the matrix AL = Ajf 8 , and the solution vector Xf = X % . However, in the 


ls The matrix Aj 3 is the transformation matrix in the forward problem in Eq. (10) initially estimated 
with the objective function /i in Eq. (17). The smoothed solution is then calculated using the function 
fs in Eq. (21) at a unique aspect ratio determined using the function / 3 in Eq. (20). The solution 
vector Xj from Eq. (17) is used to construct a trial solution as "fi = ln(W) in the case of the number 
based PSD. For the volume based PSD, the corresponding number based PSD from Eq. (21) is used 
to construct a trial solution as described in section 3 of the supplementary information. 


18 





case of the volume based PSD, the vector X+ = X? (in the case where smoothing is 
not required 19 , then the volume based PSD X\ is obtained from an objective function 
similar to f\ in Eq. (17) as described in section 3 of the supplementary information, 
otherwise Eq. (21) is used), the matrix A 7 will be scaled accordingly (see section 
3 of the supplementary information for details), and the vector C'J = C* where the 
transformed CLD C* is calculated as 

N 

C* = Y, A ^ ( 22 ) 

i=1 


where 


Xj 


* = 


Ei=i x t 


(23) 


The operator V| is a finite difference approximation to the second derivative of the 
vector Xf given as [27] 20 



h-X+ t - (h+ + h-)Xt + h+Xj ti 
\h+h-(h + + hJ) 


(24) 


where 


h_ Tj Li_i 

Li (2d) 


and X 4 + has been treated as a function of the characteristic particle size L. The pa¬ 
rameter A 2 sets the level of penalty imposed on the second derivative of Xf. If the 
value of A 2 is sufficiently large, then the penalty on the second derivative causes the LM 
algorithm to search for a solution vector which is smooth thereby avoiding solutions 
with localised oscillations (see subsection 4.4 of the supplementary information). 

The volume based PSD obtained from Eq. (21) is normalised and converted to a 
probability density distribution as 


A7 = 




(L, +1 - Li) £Zi X 


(26) 


8. Results and Discussion 

The results obtained with the two methods outlined in sections 6 and 7 are presented 
in this section. More details of the choice of parameter values are presented in the 
supplementary information. 


19 The vector Xf is defined as an exponential function of a fitting parameter similar to Eq. (18) for 
Xi. The optimisation is then performed to obtain the optimum value of the fitting parameter using 
the LM algorithm similar to the case of X 7 in Eq. (18). 

20 The form of the central difference approximation to the second derivative of the vector X+ given 
in Eq. (18) is necessary since the grid for L is non uniform as seen in Eq. (8). 
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Figure 8: (a) Experimentally measured (symbols) and calculated (solid line) CLD for PS. The calcu¬ 
lated CLD was obtained by solving the forward problem in Eq. (10) using the number based PSD 
calculated with the objective function f 3 (using X 2 = 10 -2 ) in Eq. (21). The calculation was done at 
the aspect ratio r — 1 where the objective function f 2 (in Eq. (20) using Ai ss 0.95) reaches a mini¬ 
mum. The matrix A in Eq. (10) was calculated by Method 1 (described in subsection 6.1). (b) The 
objective function f 2 (using Ai « 0.95) in Eq. (20) for different aspect ratios for PS. (c)Calculated (by 
Method 1) volume based PSD for the PS sample. The volume based PSD was calculated using the 
objective function f 2 (at Ai ~ 0.95). The objective function f 3 was not used in the calculation of the 
volume based PSD in this case as smoothing was not required. 


8.1. Results from Method 1 

Figure 8(b) shows the objective function / 2 as a function of the aspect ratio r for 
PS. The function reaches a minimum at r = 1 suggesting spherical particles. This 
is consistent with the shape of the particles in Fig. 2(a) and the mean aspect ratio 
of r « 0.8 obtained from the shape descriptor in Fig. 5(d) for this sample. This is 
also in agreement with the histogram in Fig. 6(d) which suggests that the majority of 
the particles in the sample are near spherical. Hence the aspect ratio predicted with 
Method 1 gives a reasonable description of the shape of the particles in the population 
as previously established [6]. 

The calculation in Fig. 8(b) was done with Ai ~ 0.95 in Eq. (20) (see subsections 
4.3 and 4.4 of the supplementary information for details on how the values of Ai and 
A 2 are chosen in Eqs. (20) and (21)) using the optimum transformation matrix from 
Eq. (17). Using this optimum transformation matrix and r = 1, a number based 
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PSD is calculated from Eq. (21) with the smoothness penalty set by A 2 = 10 -2 . The 
CLD corresponding to this number based PSD is shown by the solid line in Fig. 8(a). 
Furthermore, the volume based PSD (calculated at r = 1) is shown in Fig. 8(c). In this 
case, the volume based PSD was calculated from the objective function f 2 (with the 
CLD C* replaced with transformed CLD C* and the matrix A, J% rescaled as described 
in section 3 of the supplementary information) at Ai ~ 0.95. The objective function 
was not used in calculating the volume based PSD in this case as smoothing was not 
required. 

The calculated CLD in Fig. 8(a) has a near perfect match with the measured CLD 
for PS which is shown by the symbols in Fig. 8(a). The calculations in Fig 8 were done 
at a value of AT = 2 (where A f is defined in Eq. (12)). This value was sufficient to give a 
wide enough range of aspect ratios to find a good match to the experimentally measured 
CLD. If the value of A f is not large enough, then the calculated CLD may not match 
the experimentally measured CLD as the particles do not have exactly the same shape 
and Method 1 only uses a single aspect ratio to describe the shape of all the particles in 
the population. The single aspect ratio chosen will then not be representative of all the 
particles in the population. However, the imaging data narrows down the search space 
for a representative aspect ratio, and hence reduce the risk of predicting an unreliable 
aspect ratio. 

Figures 9 and 10 are similar to Fig. 8 but for COA and Glycine respectively. Figure 
9(b) shows the objective function / 2 (in Eq. (20)) with aspect ratio for COA. The 
function / 2 in Fig. 9(b) predicts an aspect ratio r = 0.3 for COA. This is reasonable 
when compared with crystals in Fig. 2(b) and the shape descriptor in Fig. 5(b). Also 
the mode of the histogram in Fig. 6(e) is close to the aspect ratio r = 0.3. The 
predicted aspect ratio of r = 0.3 in Fig. 9(b) is also close to the estimated aspect ratio 
of r « 0.2 from the shape descriptor in Fig. 5(b). Furthermore, the calculated (in a 
manner similar to the case of PS in Fig 8(a)) CLD for COA shown by the solid line in 
Fig. 9(a) has a near perfect match with the measured CLD for the sample shown by 
the symbols in Fig. 9(a). The calculations were done with J\f — 4 in Eq. (12). 

The volume based PSD for COA (calculated in a manner similar to the case of PS 
in Fig. 8(c)) is shown in Fig. 9(c). The calculated PSD in Fig. 9(c) has a left shoulder 
extending to needle lengths of about 10ywm. This gives a hint of the presence of a 
significant number of short needles in the COA sample. Some of these short needles 
can be seen in the image of Fig. 2(b). 

The calculations for Glycine in Fig. 10 are similar to the cases of PS and COA in 
Figs. 8 and 9 respectively. The predicted particle shape represented by r = 0.4 in Fig. 
10(b) is consistent with the particles in Fig. 2(c) and shape descriptor (which yields 
r « 0.4) in Fig. 5(c) as well as the histogram in Fig. 6(f). The calculated CLD (solid 
line Fig. 10(a)) also matches the experimentally measured CLD for the Glycine sample 
(symbols in Fig. 10(a)). The calculated volume based PSD for this sample at r = 0.4 
is shown in Fig. 10(c). The volume based PSD in Fig. 10(c) also has a left shoulder 
extending to about 10/jrn similar to the case of COA in Fig. 9(c). The calculations in 
Fig. 10 for Glycine were done with J\f = 2 which was sufficient to get a good match for 
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Particle length (/um) 


Figure 9: Similar to Fig. 8 but for COA. The calculation was done at the aspect ratio r = 0.3 where 
the objective function R (using Ai = 0.54) reaches a minimum in (b). The number based PSD (used 
for calculating the CLD) was calculated using A 2 = 0.05 in Eq. (21) while the volume based PSD 
(shown in (c)) was calculated using A 2 = 1CT 7 in Eq. (21). 


the measured CLD. 

8.2. Results from Method 2 

The aspect ratios in Figs. 6(b) and 6(c) show a spread over a significant range of 
particle sizes. Hence the technique referred to as Method 2 in subsections 6.2 and 7.2 
was also applied in the analysis of the data from COA and Glycine. 

The solid line in Fig. 11(a) shows the calculated CLD for COA using Method 
2. The calculation was done by searching for the optimum values of L rnm and L max 
while constructing different transformation matrices as outlined in subsection 6.2. The 
search for the optimum values of L min and L max (and hence the optimum transformation 
matrix) is done by minimising the objective function f\ in Eq. (17). Once the optimum 
transformation matrix is found, then a number based PSD is calculated by minimising 
the objective function (with A 2 = 0.1) in Eq. (21). The CLD corresponding to this 
number based PSD is shown by the solid line in Fig. 11(a). The transformed CLD C* 
in Eq. (22) is calculated using the optimum transformation matrix and the number 
based PSD obtained with Eq. (17). 

The calculated CLD in Fig. 11(a) (solid line) has a near perfect match with the 
experimentally measured CLD (shown by the symbols in Fig. 11(a)) for COA. This is 
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Figure 10: Similar to Fig. 9 but for Glycine with calculations done at r = 0.4. The values of Ai = 0.41 
and A 2 = 0.01 were used in Eqs. (20) and (21) respectively for the number and based PSD. The value 
of A 2 = 10 -6 was used in Eq. (21) for the volume based PSD. 


similar to the situation in Fig. 9(a) where the calculation was done with Method 1. 
The degree of agreement of the calculated CLD in Fig. 11(a) with the experimentally 
measured CLD demonstrates the level of accuracy that can be achieved with Method 
2. Note that the aspect ratios of each of the subgroups of each bin (in Fig. 7) were 
assigned equal weights; a simple approach that is sufficient for reasonable results in 
this case. The volume based PSD for COA (obtained by Method 2) suggests particle 
sizes from about 3/im to about 400/irn (Fig. 11(b)). This is close to the prediction of 
particle sizes from about 7/urn to about 400/zm by Method 1. Even though the ranges 
of particle sizes predicted by both methods are close, Method 2 has the advantage that 
aspect ratio is not used as a fitting parameter which removes the issue of estimating the 
optimum aspect ratio from the problem. The aspect ratio is assumed to vary according 
to imaging data available. 

Although the particle sizes estimated from these 2D images are not very accurate 
because of the focusing problem highlighted earlier, a comparison of the estimated PSD 
from the images with the volume based PSD obtained by both methods can still be 
made. This comparison shows good agreement of the estimated volume based PSD 
from the images with the volume based PSD estimated by Methods 1 and 2. Details 
are given in section 6 of the supplementary information. The peak of the volume based 
PSD obtained by Method 2 is higher than that obtained by Method 1 in this case. This 


23 














10 ° 10 1 10 2 10 3 
Chord length (/im) 



10 1 10 2 10 3 
Particle length (/im) 



10 1 10 2 10 3 
Particle length (/im) 


Figure 11: (a) Experimentally measured (symbols) and calculated (solid line) for COA. The calculated 
CLD was obtained by solving the forward problem in Eq. (10), where the matrix A was calculated by 
Method 2 as outlined in subsection 6.2. The number based PSD used in solving the forward problem 
was obtained from the objective function / 3 in Eq. (21) for A 2 = 0.1. (b) The blue diamonds are the 
calculated (by Method 1) volume based PSD for COA shown in Fig. 9(c). The black asterisks are 
the volume based PSD calculated by Method 2 for COA. The volume based PSD by Method 2 was 
calculated from Eq. (21) at A 2 = 3 x 10 -6 . (c) Similar to (a) but for Glycine. The number based 
PSD was obtained from Eq. (21) at A 2 = 0.01. (d) Similar to (b) but for Glycine, where the value of 
A 2 = 10 -5 has been used for the volume based PSD calculated by Method 2. 


is because the volume based PSD by Method 2 is slightly narrower within the size range 
of about 50yum to about 200/nn (Fig. 11(b)) so that the main peak gets higher to satisfy 
the normalisation constraint in Eq. (26). The main peak is accompanied by a smaller 
peak at a particle length close to 30/mi (Fig. 11(b)) suggesting a bimodal distribution 
for the COA particles. However, this feature of a bimodal distribution is not picked up 
by Method 1 (Fig. 11(b)). This could be because Method 2 is more efficient in picking 
out bimodal distributions in a population of particles where there is a variation of aspect 
ratio for particles of different sizes (see section 7 of the supplementary information for 
details) than Method 1. 

The situation for Glycine is similar to that of COA. The solid line in Fig. 11(c) 
shows the calculated (calculated in a manner similar to the case of COA in Fig. 11(a)) 
CLD for Glycine. The calculated CLD in Fig. 11(c) also has a near perfect match with 
the experimentally measured (symbols) CLD in Fig. 11(c). This is similar to the case 
of COA in Fig. 11(a). The calculated volume based PSD (by Method 2) in Fig. 11(d) 
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for Glycine also covers about the same range of particle sizes as in the case of Method 
1 (Fig. 11(d)) with the PSD very similar from both methods. 

The volume based PSD for COA obtained by the two Methods (Fig. 11(b)) cover 
a size range of < 10/mr to about 400/irn, while the CLD data for COA (in Fig. 11(a)) 
shows a maximum chord length of about 300/tm 21 . The PSD in Fig. 11(b) agrees 
with the image data in Fig. 6(b) for COA where the scatter plot is dense in the region 
between about 30/mr to about 300/im with a small number of particles of sizes > 300/an. 
Particles of small sizes below about 30/un are not picked up by the image processing 
algorithm because objects smaller than that are rejected by the algorithm to reduce the 
risk of processing background noise as real objects. This is the reason why particles of 
sizes < 30/trn do not contribute to the scatter plot of Fig. 6(b) even though the volume 
based PSD for COA in Fig. 11(b) suggests the presence of these particles. 

The situation with Glycine is similar to that of COA as seen in Fig. 11(d). The 
volume based PSD obtained by both methods cover a size range from about 10/un to 
about 500/un in agreement with the CLD for Glycine (in Fig. 11(c)) which shows the 
longest chord to be about 400/nn. The data in Fig. 11(d) also agrees with the image 
data in Fig. 6(c) which shows particle sizes up to about 500/un. There may be a larger 
number of large particles (of sizes close to 500/un) in the Glycine slurry than in the 
COA slurry so that their contribution to the CLD is more significant. 

9. Conclusions 

Two different methods have been developed to constrain the search space of aspect 
ratio(s) for particle size estimation using CLD and imaging data. Both methods esti¬ 
mate aspect ratio from images and then use the information in the estimation of aspect 
ratio and/or PSD from CLD data. 

In the first method, the PSD estimation from CLD data is carried out using a single 
representative aspect ratio for all the particles in the slurry. However, the search space 
for this representative aspect ratio is reduced by means of data from the images of 
the particles captured in-line during the process. This reduces the risk of predicting 
an aspect ratio which is not representative of the particles in the slurry, and hence an 
unreliable PSD. 

In the second method, a range of aspect ratios (also estimated from the images 
of the particles captured in-line) is assigned to particles of different sizes. This takes 
aspect ratio estimation out of the problem, and hence eliminates the risk of estimating 
a PSD at an aspect ratio which is not representative of the particles in the slurry. 

The techniques presented in this work have been developed to be applied in situ, 
and an in-line imaging tool has been used in this work. The currently available in¬ 
line imaging tools are not suitable (when used on their own) for obtaining accurate 
PSD and aspect ratio due to various issues outlined in the text. The limitations of 


21 Note that CLD is number based and therefore much less sensitive to presence of a small number 
of large particles. 
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using images from these in-line tools alone to get aspect ratio and/or PSD estimates 
also show up in the large error bars in Figs. 14(d) to 14(f) in subsection 4.6 of the 
supplementary information. Hence the methods presented here combine imaging and 
CLD data obtained in-line to obtain more robust estimates of PSD and aspect ratio. 
Note that the methods presented here can be applied to combine CLD with imaging 
captured with any in situ tools. The images need to be of sufficient quality so that 
aspect ratio information can be obtained from them using a suitable image processing 
algorithm. 
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Supplementary Information 


1. Probability density function (PDF) for single particle chord length dis¬ 
tribution (CLD) 

The Li and Wilkinson (LW) model gives a probability density function (PDF) which 
can be used to calculate the relative likelihood of obtaining a chord of length s from a 
particle of length L and aspect ratio r [5]. The LW model was used in this work because 
it covers the entire range of aspect ratios r G [0,1]. The PDF of the LW model is derived 
from an ellipsoid of semi major axis length a, semi minor axis length b and aspect ratio 
r = b/a [5]. For such an ellipsoid, the probability Px, ; ( 6 j,cn s j+i,a) °f obtaining a chord 
whose length lies between s t and s i+ i from a particle of length L t = 2a, depends on the 
angle a between the cutting chord and the x axis [5]. The angular dependent PDF is 
given by 
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where t = tan (a). The angle independent PDF is then given as 
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2. Determining the number of particles for making aspect ratio estimates 


Consider a situation where a total of N ob j are detected by the image processing 
algorithm. A number Nk < N ob j of objects are chosen at random from the N ob j objects 
detected, and then the average aspect ratio fkp of this sample is calculated. Then 
another random sample (with the same number of objects TV*,) is selected, and again 
the average aspect ratio fk p is calculated. This process is repeated N pic k times (at the 
same sample size) and each time the value of fk p is calculated for that sample size 
Nk- This implies that the index p takes values from 1 to N pic k. Then the sample size 
is increased up to N T , and each time the overall process is repeated. The standard 
deviation for different sample sizes can be calculated as 



Nk Nk Nk 


Figure 1: Scatter plot of the mean aspect ratio for all random selections of samples of different sizes 
for (a) PS, (b) COA and (c) Glycine. Variation of the standard deviation of the mean aspect ratio for 
all random selections of samples of various sizes for (d) PS, (d) COA and (f) Glycine. 
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The average aspect ratio (ry p ) for all random selections N pic k = 100 of samples of 
different sizes from Nk = 100 are shown in Figs. 1(a) to 1(c) for the three materials PS 












(Fig. 1(a)), CoA (Fig. 1(b)) and Glycine (Fig. 1(c)). The figures show small variations 
of the mean aspect ratio with different sample size. However, the standard deviation 
crfc consistently decreases with increasing sample size (as seen in Figs. 1(d) to 1(f)) as 
expected. The results in Fig. 1 clearly show that the calculated mean value from the 
objects detected in the images become more representative of the particles in the slurry 
as the number of detected objects increase. However, detecting more objects implies 
processing more images and the time to do this depends on the acquisition frequency 
of the image acquisition device and the image processing algorithm. Hence a decision 
needs to be made on the degree of accuracy that is sufficient for a particular process. 
Once that decision is made then the average aspect ratio can be retrieved at that sample 
size. In any case, Fig. 1(d) to 1(f) show that the standard deviation cr^ begins to level 
off at sample size Nk ~ 500. This implies that the error incurred in under sampling the 
particles in the slurry become minimal for sample sizes Nk > 500. 


3. Calculating volume based PSD 


The technique for calculating the volume based PSD is the same as that presented 
in [6] . A generalisation of the technique is necessary for the case of Method 2 (described 
in subsection 6.2 of the main text) where multiple aspect ratios are assigned to particles 
of the same characteristic length. The updated technique is described in this section. 

Obtain the normalised number based PSD X, as 


1- - X * 

* Y ’ 
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where Ah is the number based PSD calculated with the inversion algorithm. Then 
calculate the CLD C* given as 

N 

C) = V Aji X„ (8) 

i= 1 

where A 3l is the transformation matrix corresponding to the number based PSD A*. 
The CLD C* could be associated with the number based PSD X t (as in Eq. (8)) or 
the volume based PSD X? depending on the weighting applied to the matrix A 3i . The 
technique for weighting the matrix A ]t in order to associate the transformed CLD C* 
to the volume based PSD is described below. 

The volume based PSD is defined as [28] 
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where Fj is the volume of the particle with characteristic size Lj. The shape of the 
particles in this work have been approximated with ellipsoids so that the volume of 
each particle is given as 
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where a*, bi and e,; are the semi axes lengths in the x, y and z directions respectively 
of the ellipsoid of characteristic length L ls = 2a*. In this case, the origin of coordinates 
has been placed at the centre of the ellipsoid with the z direction parallel to the major 
axis of the ellipsoid. Assuming the axes lengths 6 * and q are equal, then using b{ = r^ai 
(where f, is the mean aspect ratio of all particles of the same characteristic length L t ) 
in Eq. (10) and substituting in Eq. (9) gives 
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Substituting Eq. (11) in Eq. ( 8 ) gives 
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Equation (12) is the forward problem for the volume based PSD. If the weighted 
(due to Eq. (13)) volume based PSD X t is known, then the CLD C* in Eq. ( 8 ) can be 
calculated using Eq. (12). In the case of Method 1 (described in subsection 6.1 of the 
main text) where the same aspect ratio is used for all the particles in the slurry, then 
the quantities A Jt and X i reduce to 
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Since the volume based PSD is usually not known, then an inverse problem can be 
formulated by searching for fitting parameters 7 " which minimise an objective function 
/ of the form 

M 

/ = £ 

3 = 1 

where 

X = e fi * — 1, 2, 3,..., AT. (16) 

This approach was demonstrated in [ 6 ] to correctly reproduce the volume based PSD. 
The function / is minimised using the Levenberg-Marquardth (the Matlab implemen¬ 
tation) using suitable initial trial solution for 7 ". The initial trial solution for li is 
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constructed as 7 ^ = InpQaf), where a* is the semi major axis length (defined in Eq. 
(10)) for bin i and is the number based PSD calculated from Eq. (17) of the main 

- V 

text. When a smooth solution vector X i is required, then X t is calculated from Eq. 
(21) of the main text and the corresponding X { is calculated from the same Eq. More 
details on the procedure for choosing the smoothing parameter A 2 in Eq. (21) of the 
main text will be given in subsection 4.4. 

4. Choice of algorithms parameters 

The motivation for choosing different parameter values for the algorithms used in 
this work are presented in the following subsections. 

4.1. Choice of number of particle subgroups in each bin in Method 2 

The Method 2 presented in subsection 6.2 of the main text outlines a technique 
for assigning multiple aspect ratios to particles of the same size. In this method, the 
characteristic particle size L* representing the size of particles in bin i is associated 
with multiple aspect ratios from r™ in to r™ ax (see subsection 6.2 of the main text for 
details). To achieve this, the particles in bin i of characteristic size Li are subdivided 
into N r subgroups. Each subgroup is assigned an aspect ratio 

rf G [ r ™Vr1 • ( 17 ) 



Figure 2: The mean CLD of a group of N r particles of the same size L = 100/rm but different aspect 
ratios r G [0.2, 0.7]. 

This then allows the construction of the 3D transformation matrix A from which 
the 2D transformation matrix A is constructed by averaging across slices of the 3D 
matrix as in Eq. 15 of the main text. Since the CLD for a particle of characteristic 
size L reaches a peak at a size corresponding to the width of the particle [ 6 ] , then each 
column of the average 2D matrix A (obtained using Eq. 15 of the main text) for a 
particle of size L will contain oscillations. This is due to the variations in the widths 
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Figure 3: (a) Variation of the L 2 norm in Eq. (18) for different number of particle size bins N. (b) The 
number based PSD for different number of particle size bins, (c) Corresponding CLDs for different 
number of particle size bins. 

of the particles of the same size but different aspect ratios. Hence, for more accurate 
calculations, then the value of N r needs to be chosen sufficiently large such that the 
oscillations are minimised. 

Figure 2 shows the average CLD for a group of particle of the same characteristic 
size L = 100/irn. The N r particles in each group are assigned uniformly spaced aspect 
ratios from r mm = 0.2 to r max = 0.7. The Fig. shows that the oscillations reduce as 
the value of N r increases. The oscillations become negligible at N r > 20. However, a 
value of N r = 50 was used in the calculations in the main text for more accuracy. 

4-2. Choice of number of particle size bins in Methods 1 and 2 

ft was demonstrated in [6] that the number of particle size bins N needed to get 
reasonable PSD estimate is of the order of N > 70. Since Method 1 (described in 
subsection 6.1 of the main text) uses a single representative aspect ratio to describe 
the shape of the particles (which is similar to the approach in [6]) in a slurry, then the 
number IV = 70 of particle size bins was used in Method 1. 

However, Method 2 (described in subsection 6.2 of the main text) assigns a range of 
aspect ratios to particles of the same characteristic length L. This makes it necessary 
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to search for an optimum number of particle size bins which yields accurate solutions 
and physically realistic PSDs. 

Figure 2(a) shows the behaviour of the L 2 norm defined in Eq. (19) of the main 
text (repeated here for convenience) 

lie* - C|| = s/Ti, (18) 

where 

M 

/. = £ 

3 = 1 

and 

Xi = e 7i ,i = 1,2, 3,..., AT. (20) 

The number based PSD X t in Eq. (20) is the solution vector which minimises the 
objective function f\ in (19), and the CLD C 3 is calculated from the forward problem 
described in Eq. (10) of the main text. The L 2 norm in Eq. (18) measures the degree 
of agreement between the experimentally measured CLD C* with the calculated CLD 
Cy 

The L 2 norm in Fig. 3(a) decreases as the number of particle size bins N increases, 
and begins to flatten out at about IV = 50. This indicates that the calculations become 
more accurate as the number of particle size bins increases. However, the calculated 
number based PSD X in Fig. 3(b) becomes increasingly noisy as N increases. The 
opposite situation holds for the calculated CLD C in Fig. 3(c) which becomes less 
noisy as N increases. Hence larger values of N give more accurate CLDs but the 
corresponding PSDs become more physically unrealistic. This then leads to a choice of 
N = 50 in the calculations shown in the main text. This is because the calculations 
becomes sufficiently accurate as suggested by Figs. 3(a) and 3(c). Furthermore, the 
calculations are computationally more expensive at larger values of N. 

4-3. Choice of Ai value in Method 1 

The Method 1 presented in subsection 6.1 of the main text uses a mean aspect ratio 
to represent the shape of the particles in a slurry. It involves the estimation of the 
mean aspect ratio F from images, after which the best representative aspect ratio r is 
searched from the range 

r e [f — A fa r , f + A/ay], (21) 

where ay is the standard deviation of all aspect ratios estimated from images and A f 
is the number of standard deviations chosen. Then for each aspect ratio r, a solution 
vector X (which is the number based PSD) which minimises the objective function f\ 
(given in Eq. (19)) is calculated. 

When the value of A f is made sufficiently large, then the calculated CLD C 3 from 
the number based PSD X t will match the experimentally measured CLD C* for one or 
more values of r. The degree of agreement is quantified by calculating the L 2 norm in 
Eq. (18). 
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The desired situation would have been the case where the calculated CLD Cj at 
the average aspect ratio r matches the experimentally measured CLD C*. However, 
this is usually not the case for a number of reasons. For example, the objects in the 
images are not always in focus and a number of objects are rejected because they make 
contact with the image frame. Also the particles do not all have the same shape, so 
that the estimated average aspect ratio r will not necessarily be representative of all 
the particles in the slurry. Hence it is necessary to search for a representative aspect 
ratio from a suitable range around f as given in Eq. (21). The search for a suitable 
aspect ratio is similar to the situation presented in [6] , however in this work, the search 
range has been narrowed down by means of the estimated average aspect ratio f from 
the images. 
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Figure 4: The L 2 norm in Eq. (18) for (a) PS, (b) COA and (c) Glycine. The square of the size of 
the number based PSD retrieved from the objective function f\ in Eq. (18) for (d) PS, (e) COA and 
(f) Glycine. 

Since the inversion is not necessarily carried out at the estimated average aspect 
ratio r, then a situation similar to that discussed in [6] arises, whereby the L 2 norm 
in Eq. (18) becomes nearly flat after some critical aspect ratio r*. This can be seen 
for PS (r* ~ 0.6), COA (r* ~ 0.3) and Glycine (r* a 0.4) in Figs. 4(a) to 4(c). The 
critical aspect ratio shifts to the left as the aspect ratio of the particles decreases as can 
be seen by comparing Figs 4(a) to 4(c) to Figs. 2(a) to 2(c) in the main text. 

This calls for a technique to retrieve a unique aspect ratio which is reasonable when 
compared with the shape of the particles of interest. An objective function /2 which 
imposes a penalty on the size of the calculated number based PSD X\ was introduced 
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Figure 5: The behaviour of the terms Tf and T 2 in Eq. (22) for PS at the aspect ratios indicated in 
(a) to (d). 


in [6] to pick this unique aspect ratio. The objective function f 2 is given as (same as 
Eq. (20) in the main text) 
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The term 7j is the original objective function /j in Eq. (19), and the parameter Ai 
sets the level of penalty on the PSD size which is contained in the term T 2 . The form 
of the penalty term was chosen (as discussed in [6]) because the total variation of the 
solution vector shows a general increase with aspect ratio as seen in Figs. 4(d) to 4(f) 
for PS, COA and Glycine. If the value of Ai is chosen carefully, then a unique aspect 
ratio can be found which is reasonable when compared with the shape of the particles. 
The method for estimating Ai is outlined below. 

The idea is to estimate a value of Ai such that the penalty term To just balances 
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Figure 6: The behaviour (with aspect ratio) of the critical value of Ai at which the T 2 term first crosses 
the T) term in Eq. (22) for PS. (b) The objective function f 2 in Eq. (22) with aspect ratios for PS. 
(c) The number based PSD (for PS) calculated at r = 1 where the objective function in (b) reaches a 
minimum, (d) The corresponding CLD (black line) to the number based PSD in (c). The red line is 
the calculated CLD at Ai = 0 and the symbols are the experimentally measured CLD. 


the term Tj in Eq. (22). The procedure is as follows: 

For a given aspect ratio r, obtain an initial estimate A? from 
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where T,° and T® are computed from Eq. (10) as 
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Figure 7: Same as in Fig. 6 but for COA. 

Then construct the objective function / 2 in Eq. (22) for different values of Ai from the 
range 

Ai e ..., A?(AAi) 2 , A? (AAi) 1 , A°, A^AAi, A^AAi) 2 , ..., (25) 

where the value of AAi chosen depends on the data set being analysed. The value of 
AAi = 5 was used for the PS, COA and Glycine samples in this work. Then the value 
of Ai at which the term T 2 (in Eq. (22)) just becomes equal to the term Tj (also in Eq. 
(22)) is chosen as the optimum value of Ai for that aspect ratio. The behaviour of the 
terms Tf and T 2 for different aspect ratios for PS are shown in Figs. 5(a) to 5(d). The 
situation is the same for COA and Glycine. The term T 2 just matches the term 7\ at 
a critical value of Ai (indicated as Aj) in Figs. 5(a) to 5(d). However, as Figs. 5(a) to 
5(d) show, there is a wide disparity in the values of A^ as the aspect ratio increases. 
This is seen more clearly in Fig. 6(a) which shows the variation of Aj with aspect ratio 
r for PS, and in Fig. 7(a) and 8(a) for COA and Glycine respectively. This suggests 
that the large variance in the values of A^ be removed before a meaningful average can 
be made. This is achieved by normalising the values of A^ by the standard deviation of 
the Aj; values for a particular sample and then taking the mean value as given in Eq. 
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Figure 8: Same as in Fig. 6 but for Glycine. 


( 26 ) 

^ = ( 26 ) 

q =i 

where Af r is the number of aspect ratios at which the calculations were performed and 
c?\ is the standard deviation of the A^ values. Normalising the data set (the values of 
Ai) by their standard deviation reduces the large variance in the set so that the variance 
of the data set becomes unity. The estimated Ai (by Eq. (26)) is just large enough to 
produce a shallow minimum in the objective function f 2 . Then the predicted aspect 
ratio is consistent with the actual shape of the particles. 

The estimated value of Ai ~ 0.95 from Eq. (26) for PS gives an estimated aspect 
ratio of r = 1 as seen in Fig. 6(b). This is a reasonable estimate of aspect ratio since the 
particles of PS are near spherical as seen in Fig. 2(a) of the main text. The estimated 
number based PSD at Ai ~ 0.95 and r = 1 (using the objective function f 2 in Eq. (22)) 
is shown in Fig. 6(c). The corresponding calculated CLD at Ai ~ 0.95 and r = 1 is 
shown by the black solid line in Fig. 6(d). The experimentally measured CLD for PS 
is shown by the symbols in Fig. 6(d). The calculated CLD at Ai = 0 (using f\ in Eq. 
(19)) and r = 1 for PS is shown by the red solid line in Fig. 6(d). 
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Figure 9: (a) The terms T 3 and X 4 (for the number based PSD) in Eq. (27) at various values of A 2 
for PS. (b) The number based PSD (for PS) computed from Eq. (27) at the various values of A 2 
indicated, (c) The calculated CLDs corresponding to the number based PSDs in (b). The black line is 
the experimentally measured CLD. (d) The terms T 3 and T 4 (for the volume based PSD) at different 
values of A 2 in Eq. (27). (e) The calculated (from Eq. (27)) volume based PSD for PS at the value of 
A 2 indicated. All calculations performed by Method 1 described in subsection 6.1 in the main text. 

The number based PSD for PS in Fig. 6(c) shows large particle counts at small 
particle sizes close to 10/jrn. This is because of surface roughness in the PS particles 
which contribute a significant number of short chords approximately less than 10/um to 
the CLD. The counts from these short chords are not very obvious in Fig. 6(d) due to 
the concentration of the PS sample which makes the longer chords dominate the CLD. 
However, the counts from short chords are more obvious in a more dilute system. 

Both the red and black solid lines in Fig. 6(d) have a near perfect match with the 
experimentally measured CLD (symbols in Fig. 6(d)), which shows that the estimated 
value of Ai by Eq. (26) is reasonable. 

Figure 8(a) shows that the behaviour of A^ as a function of aspect ratio, for COA 
is similar to the case of PS. The estimated value of Ai ~ 0.45 (using Eq. (26)) gives an 
aspect ratio of r = 0.3 (using f 2 in Eq. (22)) as seen in Fig. 7(b). This estimated aspect 
ratio agrees with the needle-like shape of the particles of COA as seen in Fig. 2(b) of 
the main text. The estimated number based PSD at Ai ~ 0.45 and r = 0.3 for COA 
is shown in Fig. 7(c). The corresponding calculated CLD at Ai ~ 0.45 and r = 0.3 for 
COA (shown by the black solid line in Fig. 7(d)) shows a very good agreement with the 
measured CLD (symbols in Fig. 7(d)). The calculated CLD at Ai = 0 and r = 0.3 is 
very close to the calculated CLD at Ai ~ 0.45 similar to the case of PS. This indicates 
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Figure 10: Same as in Fig. 9 but for COA. 

that the estimated value of Ai is just at the right level. 

The situation for Glycine is similar to those of the PS and COA samples. The A^ 
values also show a decrease with aspect ratio as seen in Fig. 8(a). The application of 
Eq. (26) leads to an estimated Ai ~ 0.41 and subsequently r = 0.4. The estimated 
value of r = 0.4 also agrees with the prism-like shape of the particles of Glycine. Similar 
to the cases of PS and COA, the estimated value of Ai ~ 0.41 is just at the right level. 
This is because the calculated CLD at Ai ~ 0.41 and r = 0.4 (black solid line in Fig. 
8 (d)) has a near perfect match with the calculated CLD at Ai = 0 and r = 0.4 (red 
solid line in Fig. 8(d)) and both show a near perfect match with the measured CLD 
(symbols in Fig. 8(d)). 

4-4■ Choice of\ 2 value for both Methods 1 and 2 

The number based PSDs estimated from Eq. (22) are not necessarily smooth as 
in the cases of the PS, COA and Glycine samples seen in Figs. 6(c), 7(c) and 8(c) 
respectively. They contain different degrees of oscillations. This is because the penalty 
function T 2 in Eq. (22) is not guaranteed to force the algorithm to search for a smooth 
solution. Similarly, when the volume based PSDs are estimated from Eq. (15), they 
may not necessarily be smooth. 

In order to retrieve smooth solutions a new objective function was introduced in 
section 7 of the main text. This objective function / 3 (which contains a penalty function 
that enforces smoothness) is defined in Eq. (15) of the main text but repeated here for 
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Figure 11: Same as in Fig. 9 but for Glycine. 
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where V| is the central difference approximation to the second derivative of the solution 
vector X+ as defined in Eq. (24) of the main text. In the case of a number based PSD, 
the CLD Cf = C* (the experimentally measured CLD), and A 7 = A ]t (which is the 
transformation matrix in the forward problem in Eq. 10 of the main text). For a 
number based PSD, the solution vector Xf = X t . However, in the case of the volume 
based PSD, the vector Xf = X\. In this case, the matrix AL will be scaled as described 
in Eq. (13) and CA = C* in Eq. ( 8 ). 

It is desirable to have some systematic way of estimating the A 2 parameter at which 
the term T 4 (which is the penalty term in Eq. (27)) just balances the term T 3 (which 
is the same as the function / in Eq. (15)). However (as will be demonstrated soon), 
the term P 4 does not always equal the term T 3 at reasonable values of A 2 . In most 
instances, the term T 4 just approaches T 3 . This then leads to a situation whereby the 
value of A 2 is chosen in such a way that the smooth solution does not deviate too much 
from the unsmoothed solution. In addition to this, the behaviour of the CLD Cj in Eq. 
(12) in comparison to the CLD C* in Eq. ( 8 ) is also used in the selection of A 2 . The 
C* in Eq. ( 8 ) is obtained directly from the number based PSD. If the calculation of 
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Figure 12: Similar to Fig. 9 but calculations by Method 2 (described in subsection 6.2 of the main 
text) for COA. 


the volume based PSD from Eq. (27) is correct, then the CLD Cj obtained from this 
volume based PSD using Eq. (12) should match the CLD C* from Eq. (8). Hence a 
value of A 2 is accepted if the CLD Cj calculated from the volume based PSD (obtained 
at that value of A 2 ) does not deviate significantly from the CLD C*, and the volume 
based PSD (calculated at the specified value of A 2 ) does not deviate significantly from 
the volume based PSD calculated at A 2 = 0. Similarly, the calculated number based 
PSD from Eq. (27) is rejected it its corresponding CLD (calculated from the forward 
problem in Eq. (10) of the main text) deviates significantly from the measured CLD. 

Figure 9(a) shows the behaviour of the terms T 3 and T 4 for different values of A 2 
for PS (for the case of the number based PSD) at aspect ratio 1 r — 1. The term T 4 
approaches the term T 4 for A 2 > 10 5 . The number based PSD obtained at A 2 = 10 -3 
(shown by the green asterisks in Fig. 9(b)) has about the same degree of oscillations as 
the number based PSD obtained at A 2 = 0 (red diamonds in Fig. 9(b)). At A 2 = 10~ 2 
the number based PSD (blue crosses in Fig. 9(b)) becomes smoother and has a peak 
which is close to the unsmoothed PSD obtained at A 2 = 0. However, as the value of A 2 
is increased further, the number based PSD drifts significantly from the unsmoothed 
PSD as seen by the position of the peak of the number based PSD obtained at A 2 = 10 


1 This is the case where the particles are assigned a single representative aspect ratio estimated 

using f -2 in Eq. (22). This approach was referred to as Method 1 in subsection 6.1 of the main text. 
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Figure 13: Same as in Fig. 12 but for Glycine. 

(shown by the magenta crosses in Fig. 9(b)). 

Even though the calculated CLD at A 2 = 10 (in Fig. 9(c)) still matches the experi¬ 
mentally measured CLD, the significant drift of the peak of the number based PSD at 
A 2 = 10 (seen in Fig. 9(b)) suggests that a smaller value of A 2 is appropriate to obtain 
a smooth number based PSD for this PS sample. Hence the value of A 2 = 10 ~ 2 was 
chosen for this sample. 

In the case of the volume based PSD (for PS), the term T 4 (in Fig. 9(d)) approaches 
the term T 3 at sufficiently large values of A 2 similar to the case of the number based PSD 
in Fig. 9(a). The unsmoothed (A 2 = 0) volume based PSD estimated for PS is shown 
in Fig. 9(e). Since the volume based PSD at A 2 = 0 does not contain oscillations (as 
seen in Fig. 9(e)), then the computed solution at A 2 = 0 was accepted. Furthermore, 
the CLD Cj (for PS) calculated at A 2 = 0 (red diamonds in Fig. 9(f)) matches the 
CLD Cj (solid black line in Fig. 9(f)). 

The case of COA is similar to that of PS. The term T 3 never crosses the term T 4 
within reasonable values of A 2 as seen in Fig. 10(a). The A 2 values are said to be unrea¬ 
sonable when the smoothed number based PSD begins to drift from the unsmoothed 
solution (obtained at A 2 = 0) as seen in Fig. 10(b). This is also reflected in the cal¬ 
culated CLDs which begin to deviate significantly from the experimentally measured 
CLD when the A 2 values become too large as in Fig. 10(c). 

The oscillations at the peak of the unsmoothed (A 2 = 0) number based PSD had 
been smoothed out in the number based PSD obtained at A 2 = 0.05 as in Fig. 10(b). 
The corresponding CLD at A 2 = 0.05 has a near perfect match with the experimentally 
measured CLD for COA as in Fig. 10(c). Since the peaks of the PSDs obtained at the 
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Figure 14: The shape descriptor (defined in Eq. (4) of the main text) for individual particles of 
(a) PS, (b) COA and (c) Glycine. The mean shape descriptor with error bars for (d) PS, (e) COA and 
(f) Glycine. Each error bar is one standard deviation (standard deviation of the distances from the 
centroid for all objects detected at each angular position) above and below each data point. 


higher values of A 2 (A 2 > 0.05) drift from the unsmoothed number based PSD (Fig. 
10(b)) and their corresponding CLDs deviate from the experimentally measured CLD, 
then the CLD corresponding to the number based PSD obtained at A 2 = 0.05 was 
presented in Fig. 9(b) of the main text. 

Similar to the T 4 term for the number based PSD, the T 4 term for the volume based 
PSD does not cross the T 3 term at reasonable values of A 2 (see Fig. 10(d)). The value 
of A 2 = 10 ' was chosen for smoothing the volume based PSD for COA. This is because 
Fig. 10(e) shows that the peaks of the calculated volume based PSDs drift significantly 
from the unsmoothed (A 2 = 0) volume based PSD for large values of A 2 (A 2 > 10 “'). 
Similarly, the CLD C 3 deviates significantly from the CLD C* for X 2 > 10“' as seen in 
Fig. 10(f). 

Unlike the cases of the PS and COA samples, the T 4 term for the number based 
PSD crosses the T 3 term (at A 2 ~ 0.4) for Glycine as seen 2 in Fig. 11(a). However, as 
Fig. 11(b) shows, the oscillations in the unsmoothed solution (at A 2 = 0) are already 
smoothed out for a smaller value of A 2 = 0.01 < 0.4. Similarly to the cases of the PS 
and COA samples, the peak of the smoothed number based PSD drifts from the peak 


2 This is one of the few instances where the T 4 term crosses the T 3 term in Eq. (27). 
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Figure 15: Illustration of the scheme to combine the windowing technique introduced in [6] with aspect 
ratio from images. 


of the unsmoothed PSD for large values of A 2 (as in the case of A 2 = 1 in Fig. 11(b)). 
Similarly, the calculated CLD deviates from the experimentally measured CLD at large 
values of A 2 as seen in the case of A 2 = 1 in Fig. 11(c). Hence the value of A 2 = 0.01 
was chosen for the number based PSD for Glycine. 

Also in contrast to the cases of PS and COA, the T 4 term of the volume based PSD 
is about equal to the T 3 term (for A 2 < 10~ 6 ) of the in the case of Glycine as seen 
in Fig. 11(d). However, the selection of A 2 for the volume based PSD is done in a 
manner similar to those of PS and COA. The value of A 2 = 10 ~ 6 was chosen for the 
volume based PSD for Glycine since the left shoulder seen in the volume based PSD 
at A 2 = 0 (Fig. 11(e)) becomes less obvious for A 2 > 10 -6 . Also, the CLD C 3 deviates 
significantly from the CLD C* for A 2 > 10 -6 as seen in Fig. 11(f). 

In the case 3 of Method 2, the behaviour of the terms T 3 and T 4 in Eq. (27) is similar 
to the case of Method 1. Figure 12(a) shows the behaviour of the terms T 3 and T 4 (in 
Eq. (27)) with A 2 for COA. The term T 4 does not cross the term T 3 (within reasonable 
values of A 2 ). The term T 4 becomes nearly equal to the term T 3 for A 2 > 0.1. The 
number based PSD estimated for COA at A 2 = 0.1 is shown in Fig. 12(b). It has about 
the same degree of fluctuation (at the peak) as the estimated number based PSD at 
A 2 = 0 (in Fig. 12(b)). However, as A 2 is increased to 1, the fluctuations at the peak 


3 The case where multiple aspect ratios are assigned to different particles of the same length as 

described in subsection 6.2 of the main text 
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are penalised as seen in Fig. 12(b). As the value of Xo is increased further to A 2 = 10, 
the peak of the number based PSD begins to drift from the unsmoothed number based 
PSD as seen in Fig. 12(b). At this value of Ai = 10, the calculated CLD begins to 
deviate from the experimentally measured CLD and calculated CLDs at smaller values 
of A 2 . This situation is shown in Fig. 12(c). 

The behaviour of the T 3 and T 4 terms for the volume based PSD (for COA) is 
similar to the case of the number based PSD. The T 4 term shows a near linear increase 
for A 2 < 1CT 8 as seen in Fig. 12(d). This suggests that the smoothing effect of the T 4 
term is ineffective at these small values of A 2 . However, the smoothing effect kicks in 
for A 2 > 10~ 8 . At this value of A 2 = ICC 8 , the fluctuations in the estimated volume 
based PSD only reduce slightly as seen in Fig. 12(e). However, at A 2 ~ 3 x 1CT 6 , 
the fluctuations smoothen out as seen in Fig. 12(e). As the value of A 2 is increased to 
10~ 2 , the peak of the estimated volume based PSD begins to drift from the unsmoothed 
volume based PSD as seen in Fig. 12(e). Also, the minor peak of the estimated volume 
based PSD at A 2 = 1CT 2 decreases significantly. Furthermore, the calculated CLD C 3 at 
A 2 = 10 -2 deviates significantly from the CLD C* as seen in Fig. 12(f). This indicates 
that the value of A 2 = 3 x 10 ~ 6 is more suitable for the COA sample. 

The situation for Glycine is similar to that of COA where the T 4 term does not cross 
the T 3 term as seen in Fig. 13(a). Similarly, Fig. 13(b) suggests a suitable value of A 2 
for the number based PSD to be A 2 ~ ICC 2 . Figure 13(c) also show that the calculated 
CLD for the higher value of A 2 = 1 deviates slightly from the experimentally measured 
CLD and the calculated CLDs at smaller values of A 2 . The situation for the volume 
based PSD is similar to that of COA. The smoothing effect of the T 4 term kicks in at 
A 2 > ICC 12 (Fig. 13(d)), and Fig. 13(e) suggests a suitable value of A 2 to be about 
ICC 5 . This is also supported by Fig. 13(f) which shows that the CLD C 3 deviates 
significantly from the CLD C* for A 2 > 10” 5 . 

4-5. Image processing algorithm parameters 

The tunable parameters of the image processing algorithm described in section 4 
of the main text are summarised in Table 1. The parameter names in Table 1 are the 
same as they were used in the code. The description of each parameter and its effect is 
given in Table 1. 

4-6. Shape descriptors for individual objects 

The mean aspect ratio F (used in Method 1 to construct a range of aspect ratios 
within which to search for a representative value) is obtained from the mean shape 
descriptor of the objects from images as outlined in section 4 of the main text. Although 
this mean shape descriptor is reasonably smooth as seen in Figs. 5(d), 5(e) and 5(f) 
of the main text, the shape descriptors for individual objects contain some degree of 
oscillations. Figures 14(a), 14(b) and 14(c) show the shape descriptors for the individual 
particles in the PS, COA and Glycine samples respectively. They contain different 
degrees of oscillations due to variations in the shapes of the particles. There is also a 
contribution from errors in the detection of particles’ boundaries because some particles 
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Table 1: Summary of tunable parameters of the image processing algorithm. 


Parameter name 

Description 

Default 

value 

Effect 

Step 1: Traditional s 
image. 

ipatial median filtering is i 

ised to reduce noise in the 

prefilter. median_size 

Size of the fil¬ 
ter window 

5 pixels 

The bigger the value the stronger the 
effect. Effect is similar to a low pass fil¬ 
ter but preserves sharp (high contrast) 
edges in the image. 

Step 2: A classical method of high pass filter ( 
threshold) is used to detect objects edges. 

Laplacian of Gaussian and 

ob j ect .detector. Log 
jsize 

Size of the fil¬ 
ter window 

89 pixels 

Large window sizes remove noise at the 
expense of blurred objects. 

ob j ect .detector. Log 
.power 

Gaussian 

curve stan¬ 
dard devia¬ 
tion 

0.23327 

The bigger the value, the lower the con¬ 
trast needs to be at the edges for the 
edge be detected; this does increase the 
risk of misclassifying noise as edge. 

ob j ect -detector. edge 
-threshold 

Threshold of 
detection 

251 

The higher the value, the sharper the 
edge has to be in order to be classified 
as edge. This reduces noise at the ex¬ 
pense of not detecting an edge if the 
value is too high. 


Step 3: Object processing - detected objects are filtered to discard objects 
that do not possess desirable properties. 


object 

-processing, imclose 
_size 

Size of disk 
for the op¬ 
eration of 

morphologi¬ 
cal closing 

8 pixels 

The bigger the value, the bigger the 
gaps in the object that will be filled. 
This reduces noise in the result. The 
trade-off is that the shape of the origi¬ 
nal object might be lost. 

object 

-processing.filter 
_area-minimum 

Objects area 
less than 

specified 
value are 

discarded 

900 pixels 

If an object is too small, the chance is 
that it is noise. The bigger the value, 
then the more that small objects are 
discarded. The trade-off is that poten¬ 
tially useful objects can be discarded. 

In the final stage, sc 
the coordinates of pi 

•lid objects are reduced to their boundary only, and 
ixels on the boundaries are reported. 

border 

_processing.border 

-thickness 

Desired thick¬ 
ness of the 
border 

1 pixel 

Optionally, setting this to 2 pixels will 
produce higher angular resolution out¬ 
put. However, this is often not needed. 
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are only partly in focus. There is also a contribution from impurities in the samples. 
This leads to the large error bars on the mean shape descriptors for the three samples as 
seen in Figs. 14(d), 14(e) and 14(f). The length of each error bar in Figs. 14(d), 14(e) 
and 14(f) is two times the standard deviation of the distance from the centroid. The 
size of the error bars in Figs. 14(d), 14(e) and 14(f) strongly indicates that a reliable 
estimate of the PSD cannot be made from current in-line imaging sensors (the PVM 
in this case) alone. Hence the need to combine the imaging data with CLD data as 
carried out in this work. However, the use of image data along with CLD data leads 
to significant improvements in the PSD estimates. The contribution to the size of the 
error bars by the variation in particle shape justifies the method which assigns multiple 
aspect ratios to particles of the same length rather than using an average aspect ratio 
for all particles in the slurry. 

5. Integrating windowing technique with aspect ratios from images in Method 

2 

The windowing technique developed in [ 6 ] allows the values of L min and L max to 
be estimated in-line using the bin boundaries of the CLD histogram. The quantities 
L rnin and L max (defined in Eqs. ( 8 ) and (9) of the main text) are the minimum and 
maximum sizes of the particles in a slurry. The size of a window is the number of CLD 
bins covered by that window, and the position of a window is set by the index P. For 
example, in the case shown in Fig. 15, a window of size 2 is placed at different positions 
(indicated by P = 1, 2,..., 7 in the Fig.) along the CLD bin boundaries. At the first 
window position (indicated by P = 1 in Fig. 15), the values of si and s 3 (where si 
and S 3 are the chord lengths corresponding to the boundaries covered by the window 
at P = 1) are assigned to L min and L max respectively. Similarly, at P = 2, L min = s 2 
and L max = s 4 and so on [ 6 ]. 

Consider a simple distribution of aspect ratios shown by the blue circles in Fig. 15. 
The aspect ratios in Fig. 15 were simulated using uniformly distributed random num¬ 
bers in the interval (0,1). In a real experiment, the aspect ratios will be estimated from 
images. The aspect ratios r E [r ?mn , r max \ correspond to particle sizes L E [L l ™ in , 
as indicated in Fig. 15, where and L™ ax are the minimum and maximum particle 
sizes obtained from all the images analysed. 

The chord lengths corresponding to the bin boundaries of the CLD histogram cover 
the range = 1/un to Sm +1 = 10 3 /un as indicated in Fig. 15, where M is the number 
of CLD bins. However, the particle sizes estimated from images do not cover this 
range (L l ™ in > 1/nn and L l ™ ax < 10 3 /un as illustrated in Fig. 15). This calls for 
appropriate approximations in the situations where the L min and L max values returned 
by the window do not lie entirely within the range set by L l ™ in and L l ™ ax . The following 
procedure was implemented in this work: 

1. At position P = 1 where L max < L l ™ in , assign r max to all N r particle subgroups 
(the particle subgroups were discussed in subsection 6.2 of the main text). This 
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approach was adopted because Figs. 6(b) and 6(c) of the main text suggest that 
the small particles are more rounded than the larger particles. 

2. At P = 2 where L min < L l ™ in and L l ™ in < L max < L l ™ ax , assign r max to the first 
50% of the N r subgroups and for the remaining 50%, assign aspect ratios (from 
images) corresponding to particle sizes between and L max . 

3. At P = 4 where L™ n < L min < L™ ax and L™ n < L max < L™ ax , assign aspect 
ratios (from images) corresponding to particle sizes between L min and L max to all 
N r subgroups. 

4. At P = 6 where L l ™ in < L min < L l ™ ax and L max > L l ™ ax , assign aspect ratios 
(from images) corresponding to particle sizes between L min and to the first 
50% of the N r subgroups, and assign aspect ratio f to the remaining 50% of 
the N r subgroups, where f is the mean aspect ratio estimated from images. This 
approach was adopted because Figs. 6(b) and 6(c) suggest that the larger particles 
have aspect ratio closer to the mean value than to r max . 

5. At P = 7, where L min > L l ™ ax , assign aspect ratio f to all N r subgroups. 

This approach allows the computation of the 3 D matrix defined in Eq. (15) of the 
main text. 

6. Comparison of estimated PSD from images and PSDs from Methods 1 
and 2 

Even though the objects detected from the images suffer from focusing problems as 
discussed in the main text, an estimate of the PSD from the images can be made. For 
this purpose the particles are treated as ellipsoids whose two minor axes lengths are 
equal (as described in section 3) and their semi major axes lengths are given as a. The 
number based PSD for each sample (COA and Glycine) is obtained from a histogram 
of the counts of particle lengths L = 2a using a geometrically spaced grid (with IV = 70 
bins) running from particle lengths of 1/jrn to lOOOyum. The number based PSD is then 
normalised as described in Eq. (26) of the main text. 

The number based PSDs X estimated in this way for COA and Glycine are shown 
in Fig. 16(a). The PSDs are a bit noisy due to the focusing issues in the images. The 
left tails of the PSDs terminate close to 30/iin since objects smaller than about this 
size are discarded from the image analysis as explained in subsection 4.1 of the main 
text. Hence small particles of lengths < 30/irn are not counted in this estimate of the 
number based PSD. As the small particles are not counted, then the peak of the number 
based PSD (for COA) in Fig. 16(a) is shifted to the right of the number based PSDs 
estimated for COA in Fig. 10(b) (by Method 1) and Fig. 12(b) (by Method 2). A 
similar situation holds for the number based PSD estimated for Glycine in Fig. 11(b) 
(by Method 1) and in Fig. 13(b) by Method 2. 
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Figure 16: (a) Estimated (from imaging data) number based PSD for the COA and Glycine samples 
(as indicated in the Fig.). This estimate does not include counts of particles whose lengths are less 
than about 30/im as the image processing algorithm used here is not able to distinguish between these 
small particles and background noise in the images, (b) Estimated volume based PSD for COA by 
Methods 1 and 2 (as indicated). The black solid line is the corresponding estimate from imaging data 
for the same COA sample, (c) Similar to (b) but for Glycine. 


The volume based PSD X 1 ’ is obtained from the number based PSD (before the 
normalisation in Eq. 26 of the main text) as described in Eq. (11), and then the 
normalisation in Eq. (26) of the main text is applied. The volume based PSDs for 
COA and Glycine obtained this way are shown by the solid black lines in Figs. 16(b) 
and 16(c) respectively. Similar to the case of the number based PSD in Fig. 16(a), the 
left tails of the volume based PSDs for COA (black solid line in Fig. 16(b)) and Glycine 
(black solid line in Fig. 16(c)) estimated from images are truncated at particle length 
~ 30/irn. The small fluctuations in the number based PSDs in Fig. 16(a) become 
exaggerated in the volume based PSDs in Figs. 16(b) and 16(c) due to the volume 
weighting in Eq. (11). However, the peaks of the volume based PSD estimated from 
images agree well with the peaks of the volume based PSDs estimated by both Methods 
1 and 2 for COA and Glycine as seen in Figs. 16(b) and 16(c). 
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Figure 17: Distribution of aspect ratios for the two subpopulations of particles used to test the relative 
performance of Methods 1 and 2. 

7. Performance of Methods 1 and 2 for a bimodal distribution 

The results in Fig. 11(b) of the main text suggests that Method 2 may be more 
effective for the estimation of the PSD of a population of particles with a bimodal dis¬ 
tribution of sizes and a distribution of aspect ratio for a given size. This hypothesis is 
tested here with an idealised population of particles made up of two separate subpopu¬ 
lations of different sizes and aspect ratios. The first subpopulation consists of particles 
of aspect ratio ri = 0.9 as shown in Fig. 17 and a number based distribution of particle 
lengths X tightly packed around the mean length L\ & 50/mi as shown by the black 
diamonds in Fig. 18(b). The second subpopulation consists of particles of aspect ratio 
r 2 = 0.2 shown in Fig. 17. The number based distribution of particle lengths of the 
second subpopulation is tightly packed around the mean length L 2 ~ 500/un as shown 
by the black diamonds in Fig. 18(b). The normalised (according to Eq. (26) of the 
main text) number based PSD of the population of particles is weighted towards the 
smaller particles as shown by the black diamonds in Fig. 18(b). 

Using the known number based PSD shown by the black diamonds in Fig. 18(b), 
a CLD is constructed by solving the forward problem in Eq. 10 of the main text. The 
transformation matrix A is constructed using particle sizes running from L min = 1/irn 
to L max = 1000/mi. The particles whose size fall within the range of sizes bounding the 
mean length L\ « 50/tm (shown by the black diamonds in Fig. 18(b)) are assigned the 
aspect ratio n = 0.9, while particles whose sizes fall within the range of sizes bounding 
the mean length L 2 ~ 500/mi (also shown by the black diamonds in Fig. 18(b)) are 
assigned the aspect ratio r 2 = 0.2. 

The CLD C constructed by the above procedure is then perturbed by adding a 
small perturbation C as 

C = C + C, (28) 
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Figure 18: (a) The black diamonds are the simulated CLD for the two subpopulations of particles 
whose aspect ratio distribution is shown in Fig. 17. The blue solid line is the calculated CLD for 
the two subpopulations by Method 1 at an aspect ratio of 0.9. (b) The black diamonds represent the 
bimodal number based PSD for the two subpopulations of particles whose aspect ratio distribution 
is shown in Fig. 17. The blue solid line is the calculated number based PSD (by Method 1) for this 
population of particles at an aspect ratio of 0.9. (c) The black diamonds are the transformed CLD 
(defined in Eq. (22) of the main text) for the two subpopulations of particles referred to in (a). The 
blue solid line is the equivalent transformed CLD calculated by Method 1 at an aspect ratio of 0.9. 
(d) The black diamonds are the volume based PSD corresponding to the number based PSD shown by 
the black diamonds in (b). The blue solid line is the calculated volume based PSD by Method 1 (at 
an aspect ratio of 0.9) for the population of particles whose volume based PSD is shown by the black 
diamonds. 


where 

C = eA(29) 

The parameter e = 1CD 4 controls the scale of the perturbation while the Gaussian noise 
J\f is drawn from the normal distribution with mean a = 0 and variance cr 2 = 1. 

The smallest singular value of the transformation matrix A is of order 1CT 3 in this 
case. This implies that the noise added to the CLD will be amplified by a factor of 10 3 
by the smaller singular values of the matrix A. Hence the noise scale set by e = 10“ 4 
will be amplified to about 10 _1 in the solution vector X (by the smaller singular values 
of A) which is larger than the scale of the number based PSD shown by the black 
diamonds in Fig. 18(b). Hence the noise level set by e = 10 -4 is quite significant. 

Using the simulated CLD C in Eq. (28) (black diamonds in Fig. 18(a)) and the 
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known particle size range of L mm = lpm to L max = 1000pm, then the inverse problem 
is solved by both Methods 1 and 2 described in the main text. 

7.1. Solution by Method 1 

In this case, the transformation matrix A is constructed using the known particle 
size range. The same aspect ratio is assigned to all particles of different lengths. Then 
the inverse problem is solved using the objective function defined in Eq. (17) of the 
main text. The calculation is then repeated at different aspect ratios. The use of the 
objective function in this case is to allow an estimate of the parameter Ai to be made 
for the objective function / 2 defined in Eq. (20) of the main text. Once the value of 
Ai has been estimated, then the objective function / 2 defined in Eq. (20) of the main 
text is then used to obtain a unique aspect ratio for the population of particles, and 
this procedure results in an aspect ratio F = 0.9 for the bimodal population of particles. 
The inverse problem is then solved with the objective function / 3 (defined in Eq. (21) 
of the main text) using the aspect ratio of F = 0.9 for all particles of different lengths. 
The number based PSD obtained by this procedure is shown by the solid line in Fig. 
18(b). The corresponding CLD is shown by the solid line in Fig. 18(a). 

Even though the calculated CLD by Method 1 matches the originally simulated 
CLD (Fig. 18(a)), the corresponding number based PSD does not match the originally 
simulated number based PSD in Fig. 18(b). The calculated number based PSD by 
Method 1 in Fig. 18(b) under estimates the length of the larger subpopulation of 
particles as seen in Fig. 18(b). The calculated number based PSD by Method 1 also 
underestimates the peak of the smaller subpopulation of particles. The calculated 
number based PSD by Method 1 contains a small bulge close to particle length of 2pm 
(in Fig. 18(b)) due to the noise in the calculation. 

The black diamonds in Fig. 18(c) represent the transformed CLD defined in Eq. 
(22) of the main text which corresponds to the original volume based PSD (shown by the 
black diamonds in Fig. 18(d)). This volume based PSD is calculated from the number 
based PSD by treating the particles as ellipsoids whose two minor dimensions are equal 
as described in section 3. This then allows the volume based PSD to be calculated using 
Eq. (11). The volume based PSD (shown by the solid line in Fig. 18(d)) calculated by 
Method 1 using the aspect ratio F = 0.9 and the objective function / 3 in Eq. (21) of the 
main text does not match the original volume based PSD shown by the black diamonds 
in Fig. 18(d). This is not withstanding the fact that the corresponding transformed 
CLD (shown by the solid line in Fig. 18(c)) matches the original transformed CLD 
shown by the black diamonds in Fig. 18(c). 

The mismatch of the calculated PSDs by Method 1 to the original PSDs is due 
to the fact that Method 1 uses a single aspect ratio of F = 0.9 to match the CLDs 
at the expense of predicting the wrong particle sizes. This is not surprising as there 
are no particles with intermediate aspect ratios between rq = 0.9 and r 2 = 0.2 in the 
population to enable Method 1 make a compromise. 
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Figure 19: (a) The black diamonds are the simulated CLD shown in Fig. 18(a), while the blue solid 
line is the calculated CLD by Method 2 using the aspect ratio distribution shown in Fig. 17. (b) The 
black diamonds are the number based PSD shown in Fig. 18(b), while the blue solid line is the number 
based PSD calculated by Method 2 using the aspect ratio distribution shown in Fig. 17. (c) The black 
diamonds are the transformed CLD shown in Fig. 18(c), while the blue solid line is the corresponding 
transformed CLD calculated by Method 2 using the aspect ratio distribution in Fig. 17. (d) The black 
diamonds are the volume based PSD shown in Fig. 18(d), while the blue solid line is the predicted 
volume based PSD by Method 2 using the aspect ratio distribution in Fig. 17. 


1.2. Solution by Method 2 

Since all the particles in each of the subpopulations have the same aspect ratio, 
then it implies that there is no distribution of aspect ratios for particles of the same 
size. Then the calculations by Method 2 can be simplified by collapsing the aspect ratio 
dimension of the 3D transformation matrix defined in Eq. (15) of the main text. Then 
it becomes an issue of assigning the aspect ratios rq = 0.9 and r 2 = 0.2 to the particles 
in each of the subpopulations. Although the particle size range is known in this case, 
however, the objective function f\ (in Eq. (17) of the main text) is still used so that 
an estimate of the value of the parameter A 2 (for the objective function /3 in Eq. (21) 
of the main text 4 ) can be made. Then the objective function / 3 in Eq. (21) of the 
main text is used to estimate the number based PSD and subsequently the CLD. This 
number based PSD estimated by Method 2 is shown by the solid line in Fig. 19(b) and 
the corresponding CLD is shown by the solid line in Fig. 19(a). 


4 Recall that the objective function / 2 in Eq. (20) of the main text is not needed in Method 2. 


54 























The calculated CLD by Method 2 matches the originally simulated CLD (black 
diamonds in Fig. 19(a)), while the calculated number based PSD by Method 2 correctly 
predicts the sizes of the two subpopulation as seen in Fig. 19(b). The originally 
simulated bimodal number based PSD is shown by the black diamonds in Fig. 19(b). 
The height of the main mode of the calculated number based PSD by Method 2 (solid 
line in Fig. 19(b)) is slightly lower than the height of the main mode of the originally 
simulated number based PSD (black diamonds in Fig. 19(b)) due to the effect of noise 
on the left shoulder of the calculated number based PSD by Method 2. 

Similarly, the calculated volume based PSD by Method 2 correctly predicts the 
sizes of the two subpopulations of particles as seen by the solid line in Fig. 19(d). 
The black diamonds represent the originally simulated bimodal volume based PSD. 
However, the calculated volume based PSD by Method 2 slightly under estimates the 
height of the main mode of the originally simulated volume based PSD as seen in Fig. 
19(d). The transformed CLD (in Eq. (22) of the main text) corresponding to the 
originally simulated volume based PSD is shown by the black diamonds in Fig. 19(c). 
The corresponding transformed CLD calculated by Method 2 matches this transformed 
CLD as seen by the solid line in Fig. 19(c). 



Figure 20: Left: 2D silhouette of a prismatic object (similar to the Glycine particles in Fig. 2(c) of 
the main text) in which the major (length 2a) and minor (length 26) dimensions lie in the plane of the 
paper. The object also contains another characteristic dimension shown by the line b'. The pyramidal 
cap of the object makes an angle of 9 with the horizontal. Right: The prismatic object shown on the 
left is represented by the blue solid line, while an ellipse of equivalent major and minor dimensions is 
shown by the dashed line. The chord Ci is over estimated by the ellipse while the chord C 2 is under 
estimated by the ellipse for this value of 9 = 43°. 


8. Error estimate from shape approximation 

Crystalline particles are faceted such as the case of Glycine in Fig. 2(c) of the 
main text. However, all the particles in both the COA and Glycine samples have been 
represented as ellipsoids in this work. This approximation will introduce some errors 
in the CLD calculated with the ellipsoidal model. A robust estimate of this error will 
involve a sensitivity analysis using models for faceted objects, and then investigating 
the variation of the difference between the CLDs of the faceted and ellipsoidal models 
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Figure 21: (a) The red diamonds represent the single particle CLD for a prismatic object (shown on 
the left of Fig. 20) with dimensions a = 50/tm, b = 20/rrn inclination angle 6 = 43°. The blue crosses 
represent the numerically computed single particle CLD of an ellipse of equivalent dimensions as the 
prismatic object. The black solid line is the calculated single particle CLD of the same ellipse using 
the analytic LW model defined in Eqs. (1) to (4). (b) Similar to (a) but with dimensions a = 50/irn, 
b = 10/nn inclination angle 9 = 10° for the prismatic object. 


as the shapes change. As this analysis is beyond the scope of this work, then a simple 
estimate of this error using the faceted object on the left of Fig. 20 shall be made here. 

The object shown on the left of Fig. 20 is the 2D silhouette of a 3D object with the 
shape of a prism and pyramidal caps as in the case of potassium dihydrogen phosphate 
(KDP) [29]. In this 2D representation, the major (of length 2 a on the left of Fig. 20) 
and minor (length 2b on the left of Fig. 20) dimensions lie in the plane of the paper. 
However, the prismatic object also has a third characteristic dimension indicated as 
b' on the left of Fig. 20. This third characteristic dimension is missing from the 2D 
silhouette of an ellipsoid (where its major and minor axes lie in the plane of the paper) 
of the same length (2a) and width (2 b) as the prismatic object as shown by the dashed 
line on the right of Fig. 20. This will introduce an error in the CLD estimated with 
the ellipsoidal model. The faces of the pyramidal caps are inclined at an angle 9 to the 
horizontal as seen on the left of Fig. 20. This will cause some chords such as C\ on the 
right of Fig. 20 to the over estimated while some other chord such as C 2 on the right of 
Fig. 20 will be under estimated when 9 is large enough. However, when 9 is small, the 
chord C\ will be under estimated as the ellipsoid will be completely contained within 
the prismatic object. These effects will also introduce errors in the CLD estimated 
with the ellipsoidal model. However, the size of the error will depend on the length and 
aspect ratio ( b/a ) of the particle. It will also depend on the angle of inclination 9. 

The error estimation is carried out here by constructing the CLD of both objects 
(prismatic and ellipsoid) for given values of the dimensions of the objects and inclination 
angle of the prismatic object. For a given value of a, b and 9 , a prismatic object is 
constructed and the corresponding ellipse is constructed as well. Then each object is 
rotated through a random angle € [0, 27t], and then horizontal chords are made on each 
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object at each orientation. The CLD is accumulated for 5000 rotations and then the 
CLD is normalised such that the total probability is unity. The result for the prismatic 
object for a = 50/irn, b = 20/irn, and 9 = 43° is shown 5 by the red diamonds in Fig. 
21(a). The CLD for the corresponding ellipse is shown by the blue crosses in Fig. 21(a). 
The black solid line in Fig. 21(a) is the calculated CLD (for an ellipse with the same 
dimensions as the prismatic object) from the analytical LW model defined in Eqs. (1) 
to (4). 

Of course the analytically calculated CLD for the ellipse matches that of the numer¬ 
ically computed CLD for the same ellipse as seen in Fig. 21(a). However, the ellipsoidal 
model shows some discrepancies with the prismatic model as seen in Fig. 21(a) due 
to differences in shape. The level of discrepancy is estimated by taking the L 2 norm 
of the difference between the CLDs from the ellipsoid and the prismatic object. This 
procedure gives an estimated level of discrepancy of about 6.5%. However, the level 
of discrepancy increases when a thinner prismatic object with less inclined pyramidal 
caps is used. This time the parameter values are given as a = 50/irn, b = 10/irn, and 
6 = 10°. The CLD for this prismatic object is shown by the red diamonds in Fig. 21(b). 
The numerically computed CLD for an equivalent ellipse is shown by the blue crosses 
in Fig. 21(b), while the analytically computed CLD (using the LW model) for the same 
ellipse is shown by the black solid line in Fig. 21(b). Figure 21(b) shows that the level 
of discrepancy between the CLDs form the prismatic and ellipsoidal objects increase as 
the aspect ratio and angle of inclination 9 are reduced. The level of discrepancy this 
time is about 14%. 


5 The value of 6 = 43° corresponds to KDP particles [29]. However, this value is used here to 
represent the Glycine particles in Fig. 2(c) of the main text since the Glycine and KDP particles have 
similar shapes. A smaller value of 0 = 10° is used to represent the COA particles in Fig. 2(b) of the 
main text as the pyramidal caps of the COA particles are not so pronounced as those of Glycine. The 
parameter values used in this illustration may not necessarily represent the Glycine and COA samples 
accurately, however, they show the discrepancy that could be introduced by shape differences. 


57 



References 

[1] C. Washington, Particle size analysis in pharmaceutics and other industries, Ellis 
Horwood Limited, Chichester, England, 1992. 

[2] J. Heinrich, J. Ulrich, Application of laser-backscattering instruments for in situ 
monitoring of crystallization process - a review, Chem. Eng. Technol. 35 (6) (2012) 
967-979. 

[3] P. Barrett, B. Glennon, Characterizing the metastable zone width and solubility 
curve using lasentec FBRM and PVM, Chemical Engineering Research and Design 
80 (2002) 799-805. 

[4] E. J. W. Wynn, Relationship between particle-size and chord-length distributions 
in focused beam reflectance measurement: stability of direct inversion and weight¬ 
ing, Powder Technology 133 (2003) 125-133. 

[5] M. Li, D. Wilkinson, Determination of non-spherical particle size distribution from 
chord length measurements, part 1: theoretical analysis., Chemical Engineering 
Science 60 (2005) 3251-3265. 

[6] O. S. Agimelen, P. Hamilton, I. Haley, A. Nordon, M. Vasile, J. Sefcik, A. J. Mul- 
holland, Estimation of particle size distribution and aspect ratio of non-spherical 
particles from chord length distribution, Chemical Engineering Science 123 (2015) 
629-640. 

[7] A. Ruf, J. Worlitschek, M. Mazzotti, Modeling and experimental analysis of PSD 
measurements through FBRM, Part. Part. Syst. Charact. 17 (2000) 167-179. 

[8] J. Worlitschek, T. Hocker, M. Mazzotti, Restoration of PSD from chord length 
distribution data using the method of projections unto convex sets, Part. Part. 
Syst. Charact. 22 (2005) 81-98. 

[9] H. Li, Y. Kawajiri, M. A. Grover, R. W. Rousseau, Application of an empirical 
FBRM model to estimate crystal size distributions in batch crystallization, Cryst. 
Growth Des. 14 (2014) 607-616. 

[10] N. Kail, W. Marquardt, H. Briesen, Estimation of particle size distributions from 
focused beam reflectance measurements based on an optical model, Chemical En¬ 
gineering Science 64 (2009) 984-1000. 

[11] J. Eggers, M. Kempkes, M. Mazzotti, Measurement of size and shape distributions 
of particles through image analysis, Chemical Engineering Science 63 (2008) 5513- 
5521. 

[12] M. Kempkes, T. Vetter, M. Mazzotti, Measurement of 3d particle size distributions 
by stereoscopic imaging, Chemical Engineering Science 65 (2010) 1362-1373. 


58 



[13] S. Schorsch, T. Vetter, M. Mazzotti, Measuring multidimensional particle size 
distributions during crystallization, Chemical Engineering Science 77 (2012) ISO- 
142. 

[14] S. Schorsch, D. R. Ochsenbein, T. Vetter, M. Morari, M. Mazzotti, High accuracy 
online measurements of mutidimensional particle size distributions during crystal¬ 
lization, Chemical Engineering Science 105 (2014) 155-168. 

[15] N. Kail, H. Briesen, W. Marquardt, Advanced geometrical modeling of focused 
beam reflectance measurements (FBRM), Part. Part. Syst. Charact. 24 (2007) 
184-192. 

[16] N. Kail, H. Briesen, W. Marquardt, Analysis of FBRM measurements by means 
of a 3D optical model, Powder Technology 185 (2008) 211-222. 

[17] J. Worlitschek, M. Mazzotti, Choice of the focal point position using Lasentec 
FBRM, Part. Part. Syst. Charact. 20 (2003) 12-17. 

[18] BSI, Particle size analysis - Image analysis methods, BS ISO 13322-1:2014, BSI 
Standards Publication (2014). 

[19] A. J. Izenman, Recent developments in nonparametric density estimation, Journal 
of the American Statistical Association 86 (413) (1991) 205-224. 

[20] E. F. Hobbel, R. Davies, F. W. Rennie, T. Allen, L. E. Butler, E. R. Waters, J. T. 
Smith, R. W. Sylvester, Modern methods of on-line size analysis for particulate 
streams, Part. Part. Syst. Charact. 8 (1991) 29-34. 

[21] M. Li, D. Wilkinson, K. Patchigolla, Determination of non-spherical particle size 
distribution from chord length measurements, part 2: experimental validation, 
Chemical Engineering Science 60 (2005) 4992-5003. 

[22] H. Li, M. A. Grover, Y. Kawajiri, R. W. Rousseau, Development of an empirical 
method relating crystal size distributions and FBRM measurements, Chemical 
Engineering Science 89 (2013) 142-151. 

[23] Z. Q. Yu, P. S. Chow, R. B. H. Tan, Interpretation of focused beam reflectance mea¬ 
surement (FBRM) data via simulated crystallization, Organic Process Research & 
Development 12 (2008) 646-654. 

[24] M. J. H. Simmons, P. A. Langston, A. S. Burbidge, Particle and droplet size 
analysis from chord distributions, Powder Technology 102 (1999) 75-83. 

[25] A. Vaccaro, J. Sefcik, M. Morbidclli, Modeling focused beam reflectance measure¬ 
ment and its application to sizing of particles of variable shape, Part. Part. Syst. 
Charact. 23 (2006) 360-373. 


59 



[26] T. Roths, M. Marth, J. Weese, J. Honerkamp, A generalized regularization method 
for nonlinear ill-posed problems enhanced for nonlinear regularization terms, Com¬ 
puter Physics Communications 139 (2001) 279-296. 

[27] A. E. P. Veldman, K. Rinzema, Playing with nonuniform grids, Journal of Engi¬ 
neering Mathematics 26 (1992) 119-130. 

[28] R. G. Holdich, Fundamentals of Particle Technology, Midland Information Tech¬ 
nology and Publishing, Shepshed, UK, 2002. 

[29] A. Majumder, Z. K. Nagy, Prediction and control of crystal shape distribution in 
the presence of crystal growth modifiers, Chemical Engineering Science 101 (2013) 
593-602. 


60 



